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High frequency combustion instability continues to be a major problem 
in the development and operation of rocket engines. Most mathematical 
models simulating this phenomena involve the derivation and solution of 
complex non-linear differential equations. In an effort to overcome the 
mathematical difficulties associated with the solution. of the nonlinear 
combustion instability problems, two methods of analysis .were developed. 

In investigating the problems of combustion instability in an annular 
combustion chamber, a modified Galerkin method was used to produce. a set 
of modal amplitude equations from the general non-linear partial differen- 
tial acoustic wave equation. From these modal amplitude equations, the 
two-variable perturbation method was used to develop a set of approximate 
equations of a given order of magnitude. These equations were modeled to 
show the effects of velocity sensitive combustion instabilities by 
evaluating the effects of certain parameters in the given set of equations. 
From evaluating these effects, one can ascertain which parameters . cause . 
instabilities to occur in the combustion chamber. In this analysis, it is 
assumed that in the annular combustion chamber, the liquid propellants are 
injected uniformly across the injector face, the combustion processes are 
distributed throughout the combustion chamber and that no time delay occurs 
in the combustion processes. 
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Chapter 1 


INTRODUCTION AND LITERATURE REVIEW 

During steady operation of a liquid propellant rocket engine the 
injected propellants are converted by various physical and chemical 
processes into hot burned gases which are subsequently accelerated to 
supersonic velocity by passing through a converging-diverging nozzle. The 
operation of such an engine, however, is seldom perfectly smooth. Instead 
the quantities which describe the conditions inside the combustor (i.e. 
pressure, density, temperature, etc.) are time-dependent and oscillatory. 
Such oscillations can be of either a destructive or nondestructive nature. 
Nondestructive unsteadiness is characterized by random fluctuations in the 
flow properties and includes the phenomena of turbulence and combustion 
noise. Unsteady operation of a destructive nature, on the other hand, is 
characterized by organized oscillations in which there is a definite 
correlation between the fluctuations at two different locations in the 
combustor. Such oscillations have a definite frequency and result in 
additional thermal and mechanical loads that the system must withstand. 

Unsteady operation of the destructive variety, known as combustion 
instability, was first encountered in 1940. At that time a British group 
testing a small solid-propellant rocket motor observed sudden increases 
of pressure to twice the expected level, enough to destroy a motor of 
flight weight. Since that time every major rocket development program 
has been plagued by combustion instability of some form. These 
oscillations in the combustion chamber can have several detrimental effects. 
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In some cases, particularly in solid-propellant rockets, instability 
can cause the steady-state pressure to increase to a point at which the 
rocket motor will explode. In liquid-propellant rocket chambers experi- 
encing unstable combustion, heat transfer rates to the walls considerably 
exceed the corresponding steady state heat transfer rates, resulting in 
burn-out of the walls. If the chamber can survive these effects, mechanical 
vibrations in the rocket system can cause mechanical failure or destroy the 
effectiveness of the delicate control and guidance systems. 

The phenomenon of combustion instability depends heavily upon the 
unsteady behavior of the combustion process. The organized oscillations of 
the gas within the chamber must be coupled with the combustion process in 
such a way as to form a feedback loop. In this manner part of the energy 
stored in the propellants becomes available to drive large amplitude 
oscillations. An understanding of this coupling between the combustion 
process and the wave motion is necessary in order to predict the stability 
characteristics of rocket engines. 

Combustion instability problems in liquid propellant rocket motors 
usually fall into one of three categories according to the frequency of 
oscillation. Low frequency combustion instability, also known as chugging, 
is characterized by frequencies ranging from ten to several hundred 
hertz, nearly spatially uniform properties, and coupling with the feed 
system of the rocket. This type of instability is less detrimental than 
other forms, and the means of preventing it are well understood. Low 
frequency instability will not be considered. 

A second type of combustion instability, which is less frequently 
observed, has a frequency of several hundred cycles per second. This 
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type of oscillation is associated with the appearance of entropy waves 
inside the combustion chamber. 

The third and most important form of combustion instability is 
known as high frequency or acoustic instability. As the name suggests, 
this type of instability represents the case of forced oscillations of the 
combustion chamber gases which are driven by the unsteady combustion process 
and interact with the resonance properties of the combustor geometry. The 
observed frequencies, which are as high as 10,000 cycles per second, are 
very close to those of the natural acoustic modes of a closed-ended 
chamber of the same geometry as the one experiencing unstable combustion. 
High frequency combustion instability is by far the most destructive and 
is the type to be considered by the following analysis. 

High frequency combustion instability can resemble any of the 
following acoustic modes: (1) longitudinal, (2) transverse, and (3) 

combined longitudinal-transverse modes. Longitudinal oscillations are 
usually observed in chambers whose length to diameter ratio is much greater 
than one; in this case the velocity fluctuations are parallel to the axis 
of the chamber and the disturbances depend only on one space dimension. 

For much shorter chambers the transverse mode of instability is most 
frequently observed. Transverse oscillations in rocket motors are 
characterized by a component of the velocity-perturbation which is 
perpendicular to the axis of the chamber but the disturbances can depend 
upon three space dimensions. Such oscillations can take either of two 
forms: (1) the standing form in which the nodal surfaces are stationary 

and (2) the spinning form in which the nodal surfaces rotate in either the 
clockwise or counterclockwise direction. Transverse combustion insta- 
bility, particularly that resembling the first tangential mode, has been 
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frequently encountered in modern rocket development programs and has been 
the subject of much current research. 

Historic Studies in the Problems of Combustion Instability 

Since the early 1950 T s much experimental and analytical research 
has been devoted to better understanding the phenomenon of high frequency 
combustion instability. Most of the theories presented prior to 1966 were 
restricted to circumstances in which the amplitudes of the pressure 
oscillations were infinitesimally small in the linear regime. Prominent 
among these are the pioneering studies of longitudinal instability by 
Crocco [1] as well as the studies of transverse instability by Scala [2], 
Reardon [3], and Culick [4]. A complete discussion of these theories is 
given in the work of Zinn [5] and will not be repeated here. 

Although linear theories provide the propulsion engineer with 
considerable insight into the problem, their applicability and usefulness 
in design is limited. The linear theories cannot provide answers to such 
important problems as the limiting value of the pressure amplitude 
attained by a small disturbance in the case of a linearly unstable engine, 
or the effect of a finite-amplitude disturbance upon the behavior of a 
linearly stable engine. In the latter case the result of many tests 
indicate that under certain conditions the introduction of sufficiently 
large disturbances into a linearly stable engine can trigger combustion 
instability. Another shortcoming of linear theories is the fact that 
their predictions cannot be compared directly with available experimental 
data; for, in the majority of cases, the experimental data is obtained 
under conditions in which the combustion instability is fully developed 
and in a non-linear regime. Therefore, theories accounting for these 
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nonlinearities associated with combustion instability are needed. A 
more detailed discussion of the nonlinear aspects of combustion instability 
can be found in a work by Zinn [5], 

In the field of finite amplitude (nonlinear) combustion instability, 
mathematical difficulities have precluded any exact solutions, and 
approximate methods and numerical analysis have been used almost exclusively. 
For this reason publications in this field are scarce. Notable among these 
is the work of Maslen and Moore [6] who studied the behavior of finite 
amplitude transverse waves in a circular cylinder. Their major conclusion 
was that, unlike longitudinal oscillations, transverse waves do not steepen 
to form shock waves. Maslen and Moore, however, considered only fluid 
mechanical effects; they did not consider the influences of the combustion 
process, the steady state flow, and the nozzle which are so important in 
the analysis of combustion instability problems. Nevertheless, pressure 
recordings taken from engines experiencing transverse instability reveal 
the presence of continuous pressure waves similar in form to those 
predicted by Maslen and Moore. 

One of the first nonlinear analyses to include the effects of 
the combustion process and the resulting steady state flow was performed 
by Priem and Guentert [7]. In this investigation, the problem was made 
one -dimensional by considering the behavior of tangential waves traveling 
in a narrow annular combustor of a liquid propellant rocket motor. They 
used a computer to solve numerically the resulting nonlinear equations for 
various values of the parameters involved. Due to the many assumptions 
involved in the derivation of the one-dimensional equations, the results 
of this investigation are open to question. 
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The successful use of the time-lag concept (see Crocco [1]) in the 
linear theories prompted a number of researchers to apply this model to 
the analysis of non-linear combustion instability. By considering a 
chamber with a concentrated combustion zone and a short nozzle, Sirignano 
[8] demonstrated the existance of continuous, finite-amplitude, longitudinal 
periodic waves. These solutions were shown to be unstable, however, thus 
indicating the possibility of triggering longitudinal oscillations. 

Mitchell [9] extended the work of Sirignano to include the possibility of 
discontinuous solutions. In this manner he was able to show that the final 
form of triggered longitudinal instability consisted of shock waves moving 
back and forth along the combustion chamber. Mitchell also considered the 
more realistic case of distributed combustion. 

In the analyses of Priem, Sirignano, and Mitchell the oscillations 
were dependent on only one space dimension. One of the first researchers 
to study finite-amplitude three-dimensional combustion oscillations was 
Zinn [5] whose work is an extension of the linear transverse theories and 
the analysis of Maslen and Moore. Using Crocco T s time lag model Zinn 
investigated the nonlinear behavior of transverse waves in a chamber with 
a concentrated combustion zone at the injector end and an arbitrary 
converging-diverging nozzle at the other end. In this case, it was 
necessary to extend Crocco’ s burning rate expression and transverse nozzle 
admittance relation to obtain the appropriate boundary conditions for the 
case when the flow oscillations are of finite size. As a result of this 
analysis Zinn was able to prove the existance of three dimensional 
finite-amplitude continuous waves which are periodic in time. In 
addition, he was able to prove the possibility of triggering combustion 
oscillations. An analytical criterion for the determination of the 
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stability of such waves was derived, but because of its complicated form 
and the limited capacity of available computers no specific numerical 
results were obtained. 

In more recent years other investigators such as Burstein ClO] 
have attempted to solve numerically the equations describing instabilities 
that depend on two space dimensions. Although the resulting solutions 
resemble experimentally observed combustion instability, this method 
requires excessive computer time, and studies of this type for three- 
dimensional oscillations will have to await the development of a much 
faster breed of computers. 

In a recent publication by Powell [ill, the problem of analytically 
and numerically analyzing multidimensional non-linear combustion instability 
was investigated. The problem In doing this is that a system of non- 
linear coupled partial differential equations whose solutions must 
satisfy a complicated set of boundary conditions governs the phenomena of 
combustion instability. These boundary conditions may describe the 
unsteady burning process of the wall of a solid propellant rocket motor; 
the conditions at an idealized concentrated combustion zone of a liquid- 
propellant rocket engine; or the unsteady flow of the entrance of a 
converging-diverging nozzle. Previously, in an effort to obtain analytical 
solutions to various combustion instability problems, investigators have 
been forced to simplify the original problem to such an extent that it no 
longer resembled the real problem that originally was to be solved. Powell 
proposed a method to perform a nonlinear stability analysis with relative 
ease. This method 5 applicable to both linear and non linear problems with 
complicated boundary conditions, was a modified form of the classical 
Galerkin method. The Galerkin method [11] is an approximate mathematical 
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technique which has been successfully employed in the solution of various 
engineering problems in the field of acoustics. Powell used this method 
to specifically study the non-linear behavior of combustion driven 
oscillations in cylindrical combustion chambers in which the liquid 
propellants are injected uniformly across the injector face and the 
combustion process is distributed throughout the combustion chamber. Based 
upon the results of his second and third order theories, the following 
nonlinear mechanisms were found to be important in determining the non- 
linear stability characteristics of the system: (1) the transfer of energy 

between modes, (2) the self -coupling of a mode with itself, and (3) a non- 
linear combustion mass source. Powell found that the self -coupling 
mechanism was important in the initiation of triggered instability, while 
the non-linear driving mechanism was important in the determination of the 
final amplitude of triggered instability. 

Statement of the Problem 

In this thesis, the problem of velocity-sensitive instability will 
be considered. Based upon previous work on this problem, only transverse 
oscillations will be considered due to mathematical simplicities. Also, 
the specific geometry of the combustion chamber to be analyzed will be 
annular or ring-like. The purpose of this thesis is to investigate the 
mechanisms which cause these instabilities due to the combustion process 
in a liquid propellant annular combustion chamber and attempt to state 
which mechanisms or conditions impose the greatest effect upon stability 
of combustion. 


In Chapter 2 of this thesis, the governing equations of fluid 
motion (i.e. , balance of mass and momentum) are stated. From the equations 
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the general acoustic wave equation for non-linear combustion is derived* 

In this derivation , both steady state and deviations from the steady-state 
conditions are considered and their effects incorporated into the general 
acoustic wave equation. 

In Chapter 3, the Galerkin method is used to obtain, from the 
general acoustic equation of Chapter 2, equations governing the modal 
amplitudes associated with the first two modes of transverse oscillation 
in a thin annular combustion chamber. These equations for the annular 
combustion chamber are solved numerically by the use of a Runge-Kutta 
program for various conditions. 

In Chapter 4, a set of approximate equations are derived from the 
modal amplitude equations presented in Chapter 3 by use of the two-variable 
perturbation technique. These resulting approximate equations are 
expressed both in the modal amplitude and amplitude-phase angle form. In 
this chapter, four special cases are presented for which closed-form 
solutions can be found. These four cases are (1) standing wave--no 
combustion, (2) standing wave—no gas dynamic nonlinearities, (3) 
traveling wave --no combustion, and (4) traveling wave — no gas dynamic 
nonlinearities. For problems not falling within the above categories, 
a numerical analysis is employed to solve approximate equations. 

In Chapter 5, the results contained in the previous two chapters 
are discussed and compared. Stability limits are obtained and the effect 
of neglecting various physical effects are discussed. In addition, the 
accuracy of the perturbation method is evaluated. A summary of the 
research contained in this thesis is presented in this chapter. 

In Chapter 6, a statement of conclusions is made along with 
recommendations for future research in this area. 
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DERIVATION OF THE GOVERNING ACOUSTIC WAVE EQUATION 

In order to investigate the non-linear combustion instabilities 
that occur in liquid propellant rocket engines , one must start with the 
balance laws of mass and momentum. Also, for this problem, a constitutive 
equation was formulated relating pressure and density. Mathematically, 
these principles are respectively 



+ V • (p u) = B* 


( 2 . 1 ) 


P 



-¥ 
+ u 



( 2 . 2 ) 


* *2 * 

P = a p, 


(2.3) 


where 


p - gas density 
t - time 


V - del operator of the system 


3 

3x* 


t + 




K 


A 

u - velocity of the gas 

B - fuel drop burning rate per unit volume 
p - pressure of the gas 

a - constant of proportionality (in this case - speed of 
sound ) . 
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The * representation denotes that the above physical quantities are 
dimensional. Equations (2.1) - (2.3) are based on the assumption that 
the fuel drops serve only as a source of mass for the gas phase. 
Interphase transfer of momentum and energy are neglected. 

Combining equations (2.2) and (2.3), the resulting equation is 




ft . 

-y * 

a V p. 


( 2 . 


For the physical situation depicted in Figure 1 



Figure 1. Schematic of a Liquid Propellant Combustion Chamber 
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A convenient non-dimensionalization of the variables is as follows: 


P = 


p o p (p o 


initial density of gas) 


u 


= a u 


l *■>■ 

V = ~ V 


t* = ^t 

a" 


*2 ^ 

P* = a P Q P 


* p o a 
B=4: B. 
L° 


Substituting these non-dimensional relations into equations (2,1), (2.3), 
and (2.4), the results are 


|£- + V • (pu) = B 


(2.5) 


3u , 

» at + pu 


V u 


= - $ 


( 2 . 6 ) 


p = P 


(2.7) 


where the unstarred quantities are dimensionless. 
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Dividing through by density p, equation (2.6) becomes 


3u . -*■ _ Vp 

— + u • V U - - — 
at P 


( 2 . 8 ) 


Since , 


^ = v In p , 

p 


the governing equations can be summarized as 


|£ ♦ ? . (pS> * B 


(2.9) 


+ £ ■ $ u = - $ P 
at 


( 2 . 10 ) 


P = P- 


( 2 . 11 ) 


It will now be shown that to the order of approximation inherent 
in these equations, the flow is irrotational , that is V x u = 0. To do 
this, take the curl of equation (2.10) and set it equal to zero. The 
resulting equation becomes 


$ x fe + u . $ u) = - V X v In p = 0 . 


( 2 . 12 ) 


Since the curl of any gradient is zero. This may be rewritten as 


7 x ^ x (u • V u) = 0. 

O t 


(2.13) 
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The vorticity ft is defined to be 


ft = V x u. 


(2.14) 


Thus , 


V 



at 


(V x u) 


Sft 

St 


From the vector identity 


A • V 1 = - t x (V x A) 


(2.15) 


it follows that 


-> 

u 




->-o . ->• 

V(Jgu^) - u x ft 


(2.16) 


Therefore , 


^ x (u • ^ u) = V x [V(lsu 2 ) - u x Q] . (2.17) 

Recognizing that the curl of any gradient is zero, equation (2.17) 
reduces to 

V x (u • ^ u) = - V x (u x ft) . (2.18) 

Using the vector identity 


V X (A x B) = (B • V)A - B(V • t) - (A . V)B + A(V • B) 
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equation (2.18) can be expressed as 


^ x (u • ^ u) = - [ • V) u - • u) - (u • 


+ u (V . 


(2.19) 


Therefore, equation (2.13) becomes 


M -($.$)£+ fi(V • u) + (u • V)tt - £(V ' h = 0. (2.20) 

3 t 

Equation (2.20) can now be modified by using the definition for the total 
(comoving) derivative which is 


DQ 

Dt 


3Q 

3t 


+ u 


(V ft), 


Substituting this expression into equation (2.20) and simplifying, the 
resulting equation becomes 


|| = a •($ u) - (fi V). u + u (V • ft), 


( 2 . 21 ) 


Rewriting u (V • ^) as u [V • (^ x u)] which is zero since the divergence 
of the curl of any vector is zero, equation (2.21) becomes 


|| = n .(7 u) - (fi V) • u . (2.22) 

The implications of this equation for a fluid starting from rest are as 
follows. At the initial instant of time (t = 0), the vorticity of any 
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fluid particle will be zero. Thus, the time derivative of the vorticity 
of the particle will be zero, implying that — = 0 at t = 0. Since 
ft = 0 and ^ = 0 at t = 0, it follows that ft = 0 at the next instant of 
time. By induction, it can be shown that ft = 0 for all time unless the 
velocity gradient becomes infinite for any t - 0, It is assumed in what 
follows that this does not occur and the flow is treated as irrotational. 

Since irrotationality has been proven, the velocity vector u can 
be expressed as 


u = V iJj (2.23) 

where ^ is the velocity potential. Substituting equation (2.23) into the 
left hand side of equation (2.10), the result is 

— + u ■ ^u = !^- + V(Jfu 2 ) - u x Q 
dt at 

= (V <j>) + ^ x . (2.24) 

For irrotational flow (ft = 0), the right hand side of equation (2.24) 
becomes 


v[!£ t ^ 


(2.25) 


Therefore, equation (2.10) can be written as 


" 3 ^ 

at 


+ h(^i> • v<|0 + In pj = o 


(2.26) 
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Spatially integrating equation (2.26) produces 


pt + % Vip • + In p = a(t ) 


(2.27) 


where a(t) is a function of integration. From equation (2.23), it can 
be seen that an arbitrary function of time can be added to ^ without 
affecting the result for u. Thus, a(t) could be absorbed into The 
same thing is accomplished by setting a = 0 which results in 


In p = - 


dip 

at 


- %vip ■ Vip 


(2.28) 


or 


P 


e 



+ ^ Vip • Vij)} 


(2.29) 


Thus, p and u are both known as functions of ip. From equation (2.9), the 
governing equation for ip can be written symbolically as 


+ pV 2 ip + Vip ■ % = B 

0 1 


(2. 30. a) 


P = P = 



+ % Vip 



(2 , 30 .b ) 


Rather than combining these quantities immediately, it is convenient to 
first make further simplifications based on the nature of the physical 
problem that it is desired to analyze. 
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Steady State Solution 

First, the steady state solution of equations (2.30) corresponding 
to purely axial motion will be found. Define the steady-state velocity 
potential <j> by 

ip = e$( z) (2.31) 

where z (assumed small) is the measure of the deviation of the density 
from its initial value (see equation 2.32 below). The bar notation will 
represent steady-state conditions. The steady-state burning rate w is 
defined from 


B = w(z), (2.32) 

While many other situations are possible, attention will be confined in 
the present work to the case when to = 0(e). To indicate this let 

® = ea (5 = 0(1)). (2.33) 

Thus, the burning rate B can be expressed as 

B = eg . (2.34) 


Equation (2.30.b) can now be written 


P 



2 fd$\ 2 

\dzj 


(2.35) 


Using the Taylor series expansion for the exponential function and 
retaining only the first two terms, equation (2.35) becomes 
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Substituting equations (2.31), (2.34), and (2.36) into equation (2. 30. a) 
and dividing the result by e yields 



(2.37) 


or 


dz 2 


3. P 2 (M. 

2 \dz 


d 2 <t> 

dz^ 


(2.38) 


Retaining only terms of 0(1) produces 



For simplicity, only the case of uniformly distributed combustion (i.e. 
a = constant) will be considered. Thus, integrating equation (2.39) one 
obtains 


= a z + C 
d z 1 


(2.40) 


where = u is the steady state velocity of the gas. 

At the injector (Z=0),u=0. Thus, = 0 and 


- d(f) - 

u = tt = a z . 
dz 


(2.41) 



20 


Deviations from Steady State 

It is now desired to investigate the stability of the steady 
state solution discussed above. Toward this end, an additional 
velocity potential related to perturbations from the steady state is 
defined by the equation 

ip = e [<£ + y> z, t)]. (2.42) 

A perturbation burning rate B is also defined by the equation 

B = a) + 003 « (2.43) 

It is assumed that w = 0(e) and this is indicated by defining a function 
a such that a = 0(1) and a ) = ae. Then equation (2.43) becomes 

B = e(o + to) . (2.44) 

Taking the gradient of equation (2.42), one obtains 

7\p = e[V | + V <j>] . (2.45) 


or 


Vip = e[u e z + V <j)3 . (2.46) 


From equation (2.42), the time derivative of ip can be expressed as 


dip e 3_4>_ 
~dt 3t . 


(2.47) 



Substituting the equations (2.46) and (2.47) into equation (2.30.b) 
and simplifying, one obtains 
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p = p = e 


-[e |f + he 2 (u 2 + 2u |f + V <p • V <(>)J 


(2.48) 


Expanding (2.48) in a Taylor series and neglecting terms of 0(e ) and 

higher produces the expression 


p = p = ! - e |^ + e 2 -h(n 2 + % ■ fy) -u |f + 


(2.49) 


Substituting equations (2.42), (2.44), and (2.48) into equation (2. 30. a) 
and dividing the result by e leads to 


9 ^ 


+ £ 


t . . . + 


v/li 

dt 


gi (U 2 t ^ ^ ~U ^ ^ ] 

1 - e |f t e 2 |^-Mu 2 + V4> • V<j>) - u 


3cf> 

3z 


+ V 2 <f> ] + (ue + V4» ) 

oZ / Z 


+ e 2 V(-ig(u 2 + 


*♦ • ? *>)- 5 H ♦ * (I?) 2 ^ ■ 


= a + eo . (2.50) 


Neglecting all terms of 0(e 2 ) and higher and recalling from the steady- 

state solution that u = = az and - a yields 

dz dz 
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0- 72 **=[^ 




34 - - 3 Z 4> 

3t 3z3t 


(*♦ • ? 0 ■ 


Substituting 


% |j- (.h ■ v*) 


h • ?H- 


into equation (2.51), results in 




- V^<f> + e 2 IV* * V ^ J + 2 u 


+ q 

3z3t 3t 


♦ If (*♦ - 0)] ■ - 


(2.53) 


where only terms of 0(1) and 0(e) have been retained. Equation (2,53) can 

9 2 <i> 

V'iq ■Fnn + Vi on r i irm 1 i ■Pi V\t r o m* nrr v2* = £_x + 0(e). 


be further simplified by observing that V 2 


Thus, the last term of equation (2.53) can be written 


‘ If (’ 2 * - 0 


If (0 ♦«<«> - 0) ■ • 


Since the other terms of 0(e 2 ) have already been neglected, consistency 
requires that this term be deleted and the equation be rewritten as 


3 2 d> 9 

3t^ " v * + e 


3 d) \ — 3 2 d) 

7 at r 2u + 


(2.54) 
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In this thesis, attention will be confined to transverse instability. 

For this situation 

<J> = cf>(x, y» t). (2.55) 

Therefore, equation (2.54) becomes 


3 2 <j> 

at 2 


- v 2 4> + 




ae . 


(2.56) 


To account approximately for frequency changes due to baffles, nozzle 
shapes, etc., a correction term of the form 


e 



(2.57) 


was introduced 
chosen so that 
equation for a 


into equation (2.56). This form, one of many possible, was 
the linearized form of equation (2.56) would reduce to Love s 
one-dimensional problem. This linearized form of (2.56) is 


9 2 <]) 3 2 c}> v 3 4 <t> 

St 2 " ' 3^ E 3x23tT 


0. 


(2.58) 


Thus, it can be seen that the value K will affect the acoustic frequencies. 
Physically, this is the purpose of baffles, nozzle shapes, and other 
physical parts of the combustion chamber. Therefore, inserting the 
correction term into equation (2.56), the resulting equation becomes 


3 2 <f> 

3t2 


v 2 4> + e 5 |£ + 2$<(> • v !£ 


K V' 


M ■ - 


ae . 


(2.59) 
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where K is the correction factor. This non-linear wave equation will be 
the basis for numerically and analytically investigating the transverse 
combustion stability problems occurring in liquid propellant rocket 


engines . 



Chapter 3 


DERIVATION OF WAVE EQUATIONS BASED UPON AN 
ANNULAR COMBUSTION CHAMBER 

In Chapter 2, there were no restrictions concerning the geometry 
of the combustion chamber in the derivation of the acoustic wave 
equation. In this chapter, however, a set of equations will be developed 
based upon a narrow annular combustion chamber. A typical cross-section 
for such a combustion chamber Is shown in Figure 2 below in dimensional 
and dimensionless form. 




(a) Dimensional (b) Dimensionless 

Figure 2. Dimensional and Dimensionless Form of a Circular 
Cylindrical Combustion Chamber 


In Figure 2 (a), the dimensional quantities are 

a 

r - radius of a typical point in the combustion chamber 
* 

R - inside radius of the combustion chamber 
* 

b - thickness of combustion chamber T s cross-section. 
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In Figure 2 (b), the dimensionless quantities are 

r - non-dimensional radius of a typical point 



R 


The first major assumption to be made in the geometry of the combustion 
chamber is 


| « 1 (3.1) 

which states that the circular cylinder can be thought of as a thin 
(ring-like) annulus . 

A 

Define the characteristic length L " by 

L* = R* . (3.2) 


In restricting the analysis to an annulus, a transformation to polar 
coordinates is convenient. Recall that the gradient and Laplacian 
operators in polar coordinates are 


V 


^r 9r 


+ 


e 9 9 
r 99 


t e 

z 


9 

9z 


V 2 


3 2 1 3 1 a 2 3 2 

3r 2 + r 3r + r 2 30 2 + 3z 2 


(3.3) 


The second major assumption for the simplification of the velocity 
potential is restricting 
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<t> = 4> ( 0 , t) 


(3.4) 


r = 1 . 


Therefore, using the operators of equations (3.3) on the function of 
equation (3.4), the results are 


V<f> 


-*• 3d) 

e e ae 


v 2 4> 



(3.5) 


Substituting the results of equation (3.5) into the general acoustic wave 
equation (2.58), the modified wave equation becomes 


3 2 <J> _ 3 2 <?> . _ f 3 M + 9 . a 2 < t ) v a4 ^ I - _ ae 

302 [ 3t + 1 30 303t * 3t 2 30 2 J ' (3.6) 


Now, express the velocity vector 

u = u + u (3.7) 

where u - steady-state velocity vector 

. . 

u - perturbation velocity vector . 

From the steady state solution in Chapter 2, the velocity vector was 
defined as 


= z 


di + " 

j e 

dz z . 


u 


(3.8) 
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Define the perturbation velocity vector by 


->■ i -> 

u = eV(f) 



( 3 . 9 ) 


Substituting equation (3.8) and (3.9) into equation (3.7) and using 
equation (2.23) results in 


-* 

u 



(3.10) 


To determine only the transverse velocity component of the perturbation 
velocity vector, subtract the perturbed velocity component along the 
axial (z) direction of the chamber from the total perturbation velocity 
vector. Thus , 



-> T -* 

u - u e . 
z z 


(3.11) 


In this case, since u = u(0, t) only, there is no perturbed velocity 
component in the axial direction; therefore, 

i _ 9<J> 

U t e JQ 6 0 • (3.12) 


It is now desired to find the burning rate a in terms of the parameters 
in the wave equation. To obtain this expression, assume velocity sensitive 
combustion with no history effects. Mathematically, the burning-rate 
function for velocity-sensitive combustion will be expressed by the purely 
phenomenological equation 
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o = a) n 



(3.13) 


where n is called the interaction index. 

Using the derived results for the general time-delay integral 
(discussed in Appendix A), the burning rate with history effects 
accounted for by a simple time delay is 


a = 



(3.14) 


where the subscript t represents the time delay. For simplicity, it 
will be assumed that 



Then, the burning rate can be expressed as 



where j = 0 - no time delay 

1 - time delay . 

Therefore, substituting equation (3.16) into equation (3.6), equation 


(3.6) can be rewritten 
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3 2 ij) 3 2 <j> 

3t2 3p + 


3 i}> 

at 


3^ 3 2 (fr 
3t 303t 


3 4 4» 

3t2902 


t nod 



(3.17) 


There is no closed form solution of equation (3.17) that appears likely. 
The main purpose of the present work is to determine the modifications of 
solutions of the usual acoustic wave equations that are caused by the 
presense of the nonlinear terms multiplied by e in equation (3.17). 

Thus, rather than attempt a finite difference numerical solution of 
equation (3.17), the following procedure was adopted. 

The solution is represented by the Fourier series 

4>( 0 , t) = f^(t) cos 0 + f (t) cos 26 + g^(t) sin 0 

+ g^(t) sin 26 + . . . (3.18) 


and initial conditions are chosen such that in the absence of the nonlinear 
terms, the exact solution can be formed using only the first two terms of 
the Fourier series. Because of the quadratic nature of the non-linearities, 
the second two terms in equation (3.18) represent a complete first order 
correction to the acoustic solution due to non-linear gas-dynamic and 
combustion effects. Only the first four terms in equation (3.18) are, 
therefore, retained and the approximate solution determined by this method 
is the simplest one capable of illustrating the influence of the nonlinear 



terms. The approximation can, of course, be improved by retaining 
additional terms in equation (3.18) but this is not investigated. 

Substituting equation (3.13) into equation (3.17) and using 
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the multiple angle formulas to simplify terms containing products of 
trignometric functions, one obtains 


d 2 f 


_L 


df 


dt2 ' "1 1 " dt " 2 dt ' *1 dt & 2 dt 


[‘ 


df df dg 


+ £«.+“ ZZ^ + 2e l f o ZZ^ + f i ZZ^ + So ITT 1 + g 


=1 dt J 


d 2 f 


+ Ke 1 + 2n ea) 


_ f l f 2 + g l g 2] “ 2 3 en “ [ f lx f 2x + g lx S 2x|| 


^ ^ g^ 115 dt ^ ^ +■ 1 1+*" ~ g 


d 2 g dg r df dg df dg 

X- 4- O r~ I rr 1_ J. -p i . j£- — a 2- _ + ] 

L S 2 dt 1 dt g l dt 2 dt 


d l g 


+ Ke — 2 ^ + 2 neco 


f l g 2 ” f 2 g l 


2 j ewn 


[ f lx g 2x 


It 


(d 2 f 

df o r 

dg i 

df d 

2. 

|dt 2 

+ 4 f 2 + to ^ t e|g 

1 dt 

f l dt J 


4Ke 


d 2 f r 

dt^ 


+ Jgtoe n 



I 

r 2 r 7 

i r 2 j. 2 \ i 


- *)««» g lT - f lT j 

L J 

L J 


cos 20 


jd 2 ^ 

)dt 2 


+< + 4g + to -rr-*— - e 


2 

dt 


df dg 

CT h f ] 

g l dt 1 dt 


+ 4 Ke 


d 2 g^ 

dt 2 " 


f l g l 


+ jtone 




sin 2 0 + 


cos 0 


L 


sin 0 


cone 


0. 


(3.19) 
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Equation (3.19) is a summation of terms composed of some function of 
time t and a term containing 0 variatioa Since the equation must be 
valid for all values of 0, each of the time dependent coefficients 
of the 0-terms must individually be equal to zero. Therefore, four 
ordinary differential equations governing the time-dependent modal 
amplitudes f ^ , g^, f^ and g^ emmerge from this analysis as the governing 
equations to be used for analysis of instability in an annular combustion 
chamber. These equations are 


d 2 f 

dt^ + f i 




df 



df 2 d gl dg 2 

dt 

+ 2e 

„ f 2 dt 

+ f l dt + S 2 dt + g l dt J 


d 2 f 

i- 

1 . - r 

Ke + 2neto 

_ f l f 2 + g l^ 

’2J " 2 d ena) |_ f i T f 2t + g lx g 2T 


(3. 20. a) 


d 2 g 

dti 1 + g i 


+ w 


dt 


+ 2e 


df 

g 2 at 


d §. 


+ f 


1 dt 


g 


df 

2 

1 dt 


f 



2 dt 


+ Ke 


d 2 g 

] 

dt 2 


+ 2nea) 


f l S 2 f 2 g 


J - 2 j eton j”f 


f lx g 2x 


2t S 1tJ ° 

(3.20.b) 


d 2f _ df^ 

+ 4f 2 + u ^ 


+ e 


dg 


df 


p- — l _ f — J 
dt 1 dt 


+ 4Ke 


d 

dt^ 


+ n 




eoon 




0 


(3.20.C) 
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d 2 g _ dg^ 

+ 4 §2 + “ df 


- el 


g 


df 
1 dt 


+ f 


dg^ 
1 dt 


+ 4Ke 


*-h 

dt 2 



- 


- ■■ 

ume 

5i g i. 

+ jnme 

. f lT g lT. 


(3 . 20 .d) 


In the following work only instantaneous combustion will be considered. 
Thus, the appropriate equations are equations (3.20) with j = 0. These 
equations are recapitulated below. 


d^f df 


d^ + f i + “ dt 


+ 2e 


df df 

2 — L + f — 2 

’2 dt 1 dt 


+ g 


dg i 

2 dt 


+ g 


dg r 

1 dt* 


d f 

+ Ke - v, d + 2neu) 
dt z 


f l f 2 + * 


a] 50 


(3. 21. a) 


d 2 g_ 

dt? 


dg 

1 + g i + w + 2e 


df 


dg 


g 2 dt + f l dt 


g 


df 

2 

1 dt 


- f, 


ff: 

2 dt 


d z g 


+ Ke + 2nea) f iS 2 ~ f 2 g 


[ f l g 2 " f 2 g l] 


( 3 . 21 . b) 


d 2 f _ df 

+ 4 f 2 + ^ + e 


dg df 

r L - f f 

=1 dt 1 dt 


d f 

+ 4Ke + 


[ g l 2 - f l 


- 0 


(3.21.c) 
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d 2 g 


dg 


iT +4g 2 + “dT- £ 


df dg 

r — l + f — J 

3 1 dt 1 dt 


d 2 g 

+ 4Ke ^2^" " wn£ 


S l f l 


= 0 


(3.21.d) 


The equations of (3.21) were solved numerically by the use of the 
quartic (fourth-order) Runge-Kutta method. To use this method, the 
equations of (3.21) are modified by defining the quantities 


df 


df 

2 

dt 


= a. 


dg 


dt 


L = b. 


dg^ 

2 

dt 


= b. 


(3.22) 


Substituting these expressions into equations (3.20) and solving these 
equations for the highest derivative (in this case - second order), we get 


da 

dt 


■ h - 


w(a 1 ) - 2ei|f 2 (a 1 ) + 


+ S 1 (b 2 ) ) ” 2 neaS(f 1 f 2 + g^pj /(I + Ke) 
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db 


dt 


g 1 - w( b 1 ) - 2ei^g 2 (a 1 ) + f^bg) - g 1 (a 2 ) 


- f 2^ b i^) ' 2neu)(f 1 g 2 - f^g^J / (I + Ke) 


da 

2 

dt 


’4f 2 - w(a 2 ) - ei(g 1 (b 1 ) - f 1 (a 1 ) 




g 2 - f 21 

S 1 1 


/(I + 4 Ke ) 


db , 

dt 


-4g 2 - u(b 2 ) + ei^g 1 (a 1 ) + fpbpj 


+ wne (f 1 g 1 ) 


/(l + 4Ke) 


(3.23) 


where i is the gas-dynamic index. 

By the development of a computer program incorporating the Runge- 
Kutta algorithm which can solve systems of first-order ordinary differen- 
tial equations, the eight equations (3.22) and (3.23) were numerically 
solved for the eight variables a 2 , b ±9 b 2 , f ± , f 2 , , and g 2 . 

Different cases involving varying the gas-dynamic index, interaction 
index, the correction variable ( K) , and the order term (epsilon) will be 
discussed and compared with the perturbation method of solution in a 
later chapter. In Appendix B, a sample program listing this calculation 


appears . 


Chapter 4 


TWO-VARIABLE PERTURBATION METHOD APPLIED TO THE 
ACOUSTIC WAVE EQUATIONS 

In this chapter, a set of approximate equations will be developed 
from the governing equations for the modal amplitudes (3.21), by the use 
of the two-variable perturbation method. The two-variable method is well 
suited to this type problem since one expects the solution to consist of 
sinusoidal functions with slowly varying amplitude. Applying this method, 
define two variables representing time 


5 = t 


n = et . (4.1) 

Therefore, the four modal amplitudes would now be 

f x = fjU.n) 
f 2 = f 2 U,n) 
s i = Si (5 ’ n) 

g 2 = g 2 U»n) . (4 * 2) 

By applying the chain rule of differentiation, it can be shown that 


dZ _ 9_Z 3Z 
dt 


(4.3) 


and 
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i!z * £z + 2e + e2 2?z. 

at 2 H 2 HSn 3n 2 


where Z = f^, f^ 5 g 2 respectively for each of the above equations. 

By substituting equations (4.3) and (4.4) for each modal amplitude into 
equations (3.21) and keeping terms only of 0(1) and 0(e), the resulting 
equations become 


3 2f 3f, 3f, 9 f 2 3gi 

+ f i + e[2 Jcbu + a IF + 2f 2 IF + 2f i IF + 2§2 IF 


3g 9 3 z f, 

+ 2g l IF + K W + 2n “ (f l f 2 + g l g 2 ):i= 0 


3 2 g, 3 2 g -i _3g, 3f, 3g 2 df 2 

IF~ + S 1 + e[2 3?3n + °W + 2§ 2 JF + 2fl IF ~ 2Sl IF 


3 g*i 3 2 ga _ 

- 2f 2 IF" + K "IF" + 2n “ (f l g 2 - f 2 S l^ = ° 


3 2f 2 3 2 f 2 _ 3f 2 3g x 3f x 

3? 2 + 4f 2 + e[2 a£3n + G IF + § 1 IF " f l35~ 


g2f 

+ 4K + h wn(g 1 2 -f 1 2 )] = 0 


9 2 S 2 ^2 ^1 ^^1 

3 F“ + 4§ 2 + e ^ 2 3^3n + a W ~ §1 W ~ f;L HT 


+ 4Kj jr- - ujn(f lgl )] = 0 . 
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From the straight-forward perturbation method, define the modal amplitudes 
by the series expansions 

f- L = f 1Q (£>n) + e f i:L (5,n) + • 

f 2 = f 20 (5,n) + e f 21 (?,n) + • 

g i = g io (?,n) + e g n ( S’ n) + * 

s 2 ~ g 2 0 ( '^’ r,) + E g 21^’ n) + ■ 

Again by applying the rules of differentiation, it 

3Z _ 3T + 3K 

h ~ n e h 

3 2 Z 3 2 T 3 2 K 

3 ? 2 “ £ 

£2- - 2!l_ + e l!*- (4.7) 

3£3n 9£9n 9£9o 

where Z = f , , f 2 , g^> g 2 

T ~ f 10’ f 20’ g 10’ g 20 

and K = f ±1 , f 21 , g^, g 21 , respectively. 

Substituting the expressions of (4.6) and (4.7) into equations (4.5) and 
keeping terms only of 0(1) and 0(e), the resulting equations become 


. . . (4.6) 

can be shown that 
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3 2 f 


35 


10 — 9f 10 

7^ + f in + zlTTT^- + + 2 + * ^ + 2f 


3f 


10 


'10 u 3 ^ 


11 35 3 n 35 


20 35 


+ 2f 


10 3 ? 


*^20 3g 10 

+ 2gon TF~ + 2g 


3 20 35 


10 35 


3g 20 „ 32f 10 
+ K- ~ZTZ 


35' 


+ 2nw(f 1() f 20 + g 10 g 20 )l 0 


3 2 g 


10 


3 2 g 


11 


35 


7 + S 10 + e ^35 7 + S 11 + 2 353n 


32g 10 . - 3g 10 


3f 


+ a 


35 


+ 2g 


10 


20 35 


. 2f ff20 2 ^20 _ IflO 3 S 10 

+ 2r, „ « »• 2g n n 2f on T5 7 


10 35 


3 10 35 


'20 35 


+ 2nu)(f 10 g 20 f 20 g 10^ ° 


3 2 f„ 


3 z f , 


2jt [i-P 
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3f 


3 2 f, 
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- f._ + 4K -ss^ - am(f 10 g 10 )] = 0 . 


10 35 10 35 


35' 


(4.8) 


By separating the terms of 0(1) and 0(e) in the equations of (4.8) and 
equating both sets of terms equal to zero, the resulting equations become 


3 2 f 

35 1 


10 


+ f 


10 


0 



40 


a 2 g 


35 


T 2 - * *10 = 0 


a 2 f 


20 


35‘ 


+ , * f 20 = 0 


a 2 g 


35 


^ ■ 0 


s2 ^ii + f „ « . -fi» . 2f j!i° . 2f , 


35 


'11 ^353n “ a 35 


' 20 a^ 10 35 


^®10 9g 20 ^ 2 ^10 — 

~ 2g 20dl 2g 10 Ti 2 K 7I 2 2n ai( f iQ f 20 + g lo S 20^ 


3 2 g 


11 


3 2 glO _ 3S 10 


3f 10 „„ 9S 20 


352 + g ll = “ 2 35 3ri " 0 35 " 2g 20 35 " 2f 10 35 


3^20 9g 10 S2g 10 — 

+2g 10 ?5 + 2f 20 35 K 752 2n a)( f io g 2o " f 20 g 10^ 


3 2 f 


21 


32f 2o - 32f 20 3 %0 . „ 3f 10 


Ji 2 + 4f 21 ~ " 2 353n a 35 " g 10 35 + f 10 35 


3 

■ 4K ' ** ( « m 2 ■ f io 2) 


3 2 g 


3 2 g. 


3g. 


35 


21 _ o 20 - & 20 

1 — + 4 Soi - “ 2 T7F1 a — — + g 


3f 


10 


353n 35 


10 35 


3g 


+ f 


10 


3 2 g 


20 . 


10 35 4K 35 2 "" + nw(f 10 g 10 ) • 


(4. 9. a) 


(4.9.b) 



The equations of (4. 9. a) are linear second-order differential equations. 
Therefore, it can be shown that assuming the appropriate form of a solu- 
tion, the results become 
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f 10 = A 1 (n) cos 5 + B 1 Cn) sin £ 

S 10 = A 2^ COS ^ + B 2 ( ' n ^ sin ^ 
f 2 o = Ag(.n) COS 7 \ + B 3 (n) sin 2£ 

S 20 = A 4^ h) cos 2£ + B 4 Cn) sin 2£ . (4.10) 

Substituting (4.10) into (4.9.b) and using the multiple-angle formulas 
yields 
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+ 2K[-4A 4 ] - -|nu)[-|<A 1 A 2 - B^X] 


cos 2£ t 


(4.11) 


where + . . . indicates terras multiplied by sines and cosines of integral 
multiples of £ other than those shown. The particular solutions corre- 
sponding to the terms shown on the right-hand sides of (4.11) will contain 
terms proportional to £ sin n£ or 5 cos n£ tn = 1 for (4. 11. a, b), n - 2 
for (4.11.C, d)3. Thus, the second approximation would be unbounded for 
large £ while the first approximation is bounded for all 5. These 
unbounded terms are called singular terms. The terms on the right-hand 
sides of (4.11) indicated by + . . .do not lead to singular terms. 

The idea of a perturbation solution is that higher order terms in 
the series solution represent small corrections to this first term to 
obtain a uniformly valid expansion. The presence of this singularity 
causes this fundamental idea to be violated. Therefore, since the expres- 
sions of q dependency are independent of the variable causing the singu- 
larity, the n-dependent expressions can be set individually equal to zero 
to avoid this problem. Therefore, from equations (4.11), the resulting 
equations, which are eight ordinary first-order differential equations 
having t \ dependency , become 
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— + 2° A 1 + 2*®! + 2*- A l A 3 + B 1 B 3 + A 2 A 4 + 


+ -inwCB^Ag - + B 2 A^ ^2 B 4"^ ® 

dB -i -i i n _ 

dT + ? B 1 - X + X B 3 - B 1 A 3 - A 4 B 2 + SV 


V L v]mi^ + 


4 - 

4 - 


pq 


o 


PQ 

CM 

PQ 

I 

co 

< 

CM 

< 

\ 


O 

II 


CO 


CO 

< 

CM 

< 


o 

II 


4 - 

pq 

CM 

PQ 


M 


PQ 

pq* 


CM 


rH I CM 


CM 

< 

Hd 


CM 

pq 

CO 

< 

I 


pq 


4* 

zt 


pq 


< 

CO 

4 * 

zt 

( — 1 

pq 

< 

< 

< 

CM 

hTcm 

CM 

< 


< 

4 * 

+■ 

4- 

f — 1 

4 - 

CO 

CM 

pq 

CM 

pq 

i — | 

pq 

¥ 

4 * 

< 

< 

pq 

rH |CM 

4- 

H |CM 

H- 

-h 

4 " 

1 


p 

T3 


pq 

rH 

< 

I 

fa 

nfcM 

4- 


CM 


rH CM 


CM 

pq C 

Hd nd 


co 

pq 

CM 

pq 

i 

CM 

< 

CO 

< 

I 

4 - 

pq 

i — I 

pq 


< 

ti 

rHfcM 

H- 



4- 

1 — l 


rH 

CM 

pq 

CM 

rH 

pq 

< 


I 1 


CM 


CM 

CM 


PQ 

< 


CM 

1 l 


< 

1 

8 

o 

hT* 


ii 


4 - 

rn 

H- 

co 

CM 

co 

PQ 

PQ 

< 

S 

CM 


4 - 

< 

4 ~ 

4 * 

1 

1 

CO 

H 

CO 

.< 

PQ 


1 P 

i— i 

1 p 

rH |CM 

< 

rH | CM 


i i 


4 - 

i§ 

4- 

CO I 

rH fcO 

cO I 

< 1 P 


PQ P 

d |nd 

+ 

HZ) |Hd 


o 

i — i 
CM 

It 

< 

rH 

< 

1 — 1 

CM 

1 

i — 1 

pq 

CM 

4- 

pq 

rH 

CM 

rH 

pq 

hY* 

< 

• 

4* 

4* 

CM 

pq 

CM 


pq 

zf 

i 

H- 


CM 

4* 

CN 


< 

P 

1—1 

rH |CN 

1 ? 

-h 

HUD 


M 


4* 

HZ) |h3 



I — 

i 



rH 


pq 

CM 


< 



4 

i — 1 


< 

CM 


pq 

O 

u-i 


H i 

4" 

II 

1 


ro 

• — I 


4* 

PQ 

<3 

G 

CM 


/ 

< 

4- 

4* 

1 

i 

CM 


4- 

pq 

,pq 

rH 

1 

? 

< 
1 1 

H| 

[CM 


4- 

rHfoO 

4* 



PQ 

P 

4- 

nd 1 

Hd 


45 


- ^ a i a 2 - W = 0 • ^ 

Since equations (4.12) are first-order nonlinear ordinary differential 
equations, the fourth-order Runge-Kutta program, previously developed, 
can be used to solve for the modal amplitude coefficients. By finding 
these coefficients for various points in time, a relation between the 
results of equation (3.21) and equation (.4.12) can be observed to the 
approximation of order e. 

Solving equations (4.12) for the highest derivative (first order 
in this case) and substituting n = et, the governing equations for the 
Runge-Kutta program become 

-i -i 

dT =£ [ -^ A l - ^ B 1 " 2 (A 1 A 3 + B 1 B 3 + A 2 A 4 + B 2 V 

- - A x 3 3 + B 2 A 4 - A 2 B 4 )] 

dB t -i i i N 

ST = £ ♦ ^ A ! - 2< A 1 B 3 - B 1 A 3 “ \ B 2 * W 

- —nija(A^A^ + + ^2^4 + ^2^4^ 

dA 0 nil n 

JC = £ c- T k 2 - |<B 2 - 2 (a 1 a 4 + b iB4 - a 2 a 3 - b 2 b 3 ) 

- ^nwCA,^ - A x B 4 + A 2 B 3 - A 3 B 2 )3 

= e C- ^aB 2 + 7 KA 2 ~ T ( ‘ B 4 A 1 “ \ B 1 + A 3 B 2 " A 2 B 3 ^ 
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- iumUjA,, t BjB 4 - A 3 S 2 - 


dA 

dt 


ZT- = U- 1° A 3 - 4KB 3 - ¥ h 2 2 - ®2 2 + B l 2 ' A l 2) 


^“ CA i B l - A 2 B 2 )] 


dB 

dt 


- = e[- 4 iB + 4KA 3 - ^AjBj - AjB^ 


- ll"“ (A 2 2 ' ® 2 2 - A l 2 + B l 2)3 


dA l 

= e[- - 4KB 4 " 4 (B l B 2 ' A lV 


2 4 


- ■| r> w ( A 1 B 2 + A 2 B 1 ^ 


dB L 

dt 


“ = e[- |oB + 4 KA 4 + ^B 2 A x + A 2 B 1 ) 


+ g Wn ^ A ]_ A 2 ~ B 1 B 2 ^ ^ 


(4.13) 


It is often convenient to express the equations for and B^ 
terms of amplitudes, C., and phase angles, which are also functions 

of the slow time variable n. Mathematically, we can express the relation- 
ships between the quantities as 

(4. 14. a) 

A. = C. cos A. 
l i i 


B. = C i sin <|> i 


( 4 . 14 .b) 


dA. 

x 

dn 


cos 
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dC. 

1 

dn 


d$. 

d>. - C. - — sin <J>. 
x 1 dn i 


(_4 . 14 . c ) 


dB. 

i 

dn 


dC. 

l 

dn 


sm 


deb . 

) . + C . -3 cos 9 . 

i i dn i 


(4.14.d) 


where i = 1, 2, 3, and 4 for each of the equations above. Substituting 
the expressions of (4.14) into the first two equations of (4.12), the 
resulting equations become 


dC <4 ■)_ _ ! 

cos *1 - C 1 d^r sin *1 + 2 ^1 cos +1 + 2 *° 1 S1T1 *1 

t ■|cc i c 3 cos 4> 1 cos <j> 3 + CjCg sin ^ sin $ 3 


+C 2 C 4 cos <J» 2 cos «(» 4 + C 2 C 4 sin 0 2 sin <)> 4 ] + 



sin <j> cos <{> 3 - C^C 3 cos <j>^ sin (f 3 + C 2 C 4 cos (j) 4 sin <j) 2 


- C 2 C 4 cos <j> 2 sin <j> 4 ] = 0 

dC. d(|) ]_ q 

— sin 4*1 + c i dTT - cos *i + 2° c i sin " 2*^1 cos *1 

+ ^CC 1 C 3 cos $ 1 sin <p 3 - C 1 C 3 cos <|> 3 sin ^ 

- C 2 C 4 cos <f> 4 sin <), 2 + C 2 C 4 cos <J> 2 sin <f> 4 ] + -gnuItC^Cg 

cos <j) 1 cos 4> 3 + sin ^ ^ sin <ji 3 + C 2 C 4 cos cos $ 4 
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t C 2 C 4 sin $ 2 sin ^ = 0 » (4.15) 

Multiplying the first equation by cos <f>^ and the second equation by sin<J >^ 9 
adding the two expressions together, and using appropriate multiple-angle 
identities from trigonometry, the resulting equation for becomes 

dC l 1 - 1 

— + 2° C 1 + 2< C l C 3 Cc0s(2 h " V 3 
+ C 2 C 4 [cos((f) 2 - <}> 4 + 0 1 )3} + 

sin( 2<J > 1 - <{> 3 ) + C 2 C 4 sin (* 2 - <t > 4 + 4' 1 )> = 0 . (4.16) 

Similarly , multiplying the first equation of (4.15) by -sin ^ and the 
second equation by cos adding the two expressions together, and using 

appropriate multiple-angle identities for trigonometry, the resulting 
equation for <|> becomes 


-| 9 4 

— ■ 2* • 7*3 sin(2 h - V + — sinU 2-Vh )] 

c c 

+ ^nw[C 3 cos( 2(|) 1 -4> 3 ) + 2 -~ cos(<j> 1 +4» 2 — 4» 4 >] = 0 • (4.17 

Using these procedures discussed above, equations for $ 2 5 ^ 3 5 $35 ^ 4 * 

and ((> can be derived. Thus, these transformed equations are 

dC 1 

+ 2° C 2 + 2 CC 1 C 4 °°s( + i-* 4 +*2 ) - C 2 C 3 cos(2* 2 -+ 3 U 
+ -jnwCC^C^sinC <J ) 4 +4 > 2) - C 2 C 3 sin( 2<|) 2 -<|> 3 )] - 0 
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^2 1 1 ^1^4 

— " 2 K ‘ sin( ^VV - C 3 sin ( 2 VV J 


1 - C 1 C 4 

+ -^ncdE— — cos(cfi 1 -<j> 4 +4> 2 ) + C 3 cos(2<() 2 -4> 3 )I] = 0 


dC 


3 + icC 0 + ■^{C 0 2 cos( 2ef> 2 — c}> 3 ) - C J 2 cos(2({. 1 -<t.^)J 


dn 2 3 8 2 


r l r 3' 


1 njj[C 0 2 sin( 2<j> 2 -^ 3 ) - C^sinC 2^-^)] = 0 


16 u 2 


1 Y 3 J 


d h 1 C 2 C 1 

dT " 4K + 8^07 sin < 2 VV - c r 8in < 2 *i-+3 } 3 


+ I t^cT" c ° 3(24* 2 ^ 3 ) - cT" cos(2<j> 1 -<j) 3 )] = 0 
o o 


dC 


d^ + ^° C 4 " -^tC 1 c 2 cos(4> 1 +4> 2 -4> 4 )3 + ^“^2 


sin(<() 1 + <J> 2 - (}> 4 )] = 0 


^4 1 <3 1^'2 1 - ^ 1^2 

_ 4K + -if sin(4i +d> — d> )] - — n^f— — - 

dfi 4 1 - C u ^V J o 1 L C 

4 H 


cos^ + <{> 2 - <f> 4 )] = 0 . 


(4.18) 


Equations (4.16), (4.17), and (4.18) are the general combustion 
equations in terms of amplitudes and phase angles. From this point, 
special cases can be investigated isolating certain conditions and closed- 
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form solutions can be obtained for these cases. It is convenient to do 
this in order to check the closed-form results of the special cases with 
the results from the general equations (4.16), (4.17), and (4.18) when 

the same conditions are imposed. 

The first case to be evaluated is the case for standing waves 
with no combustion effects. To simulate standing wave effect, set the 
amplitudes C 2 and C 4 and phase angles <j> 2 and <j> 4 equal to zero. This 
automatically satisfies four of the eight equations (4.18). To achieve 
the no-combustion effect, set the interaction index, n, equal to zero. 
Also, set the correction variable, K, equal to zero since the effect of 
K will be investigated separately at a later time. Imposing these con- 
ditions, the governing equations reduce to 

+ I ~ aC l + l C 1 C 3 cos(2* 1 -* 3 ) = 0 (4.19.; 


d< h 1 

- 2 K 3 sin(2 ( , 1 -* 3 ): = 0 
dC o -i i 

JT * T° C 3 - I C 1 2oos(2 +1 -* 3 ) = 0 

d< K 1 C 1 2 

— ■ s = °- 

The initial conditions imposed for this case are 


( 4. 19 .b ) 


(4.19.c) 


(4.19.d) 


C^O) = 1 
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Cg(0) = 0 
<p 1 ( °) = 4> 10 
h (0) = ^30 . 


To attempt a closed-form solution, let 


C = e ' 2CrT1 F 
1 1 


(‘ 


C = e"^ n F 
u 3 3 


dC 


1 _ --%an 


dn 


= e 


V" dF l 

(- 5>r, + 


dC, i — dF 

3T ■ e ' %0 " ( - “ ,r 3 + *" % ^ar> • 


Substituting these expressions into equations (4. 19. a) and (4.19.c) 
dividing through by e , the resulting equations become 


| 003(2^ - (s )=-^ r i r 3 ■ 0 


( 



- 003(24, 



( 


Multiplying equation (4.2 2.a) by 1/4 and equation (4.22.b) by Fg/F^. 
adding the two equations, terms containing the cos( 24 1 - 4 g)e 
eliminated. In doing so, the result becomes 


(4.20) 

.21. a) 

.21. b) 

4.21. c) 

4 . 21. d) 

and 

4. 22. a) 

4. 22 . b) 
and 


are 
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fl 3 df 3 _ n 

do + 4 r x dn 0 


(4.23) 


Multiplying through equation (4.23) by F gives 


1 5T [F i 2 + “V- 1 = 0 • 


(4.24) 


Integrating with respect to n then dividing by 1/2, the resulting equa- 
tion becomes 


F 2 + 4F 3 2 = D x (4.25) 

where is a constant of integration. This constant depends upon the 
initial condition imposed on the problem. From the initial conditions 
given in (4.20) and using the transformation (4. 21. a) and (4.21.b), it can 
be shown that F^(0) = 1 and F^(0) = 0. Therefore, equals to 1. Thus, 
equation (4,25) becomes 

F ^ - i _ 4F 2 (4.26) 

o • 

Taking equation (4.26) and substituting into equation (4.22.b), then 
separating variables, the resulting equation becomes 

dF !- 

= -re ^ ffn cos(2<j> )dr) (4.27) 

0-4F 3 2] 8 1 3 

Letting 2<$> - <p = £tt, which satisfies equations (4.19.b, d), yields 

-L O 

t 

cos (2^)^ - <f> 3 ) = (-1) where l = 0, 1,2,3. . . Substituting this expres- 
sion and integrating the above equation, the resulting equation becomes 



? tanh' 1 2F 3 .it- | e^opC-l/ 


(4.28) 


where is a constant of integration. Using the initial condition F^(0)- 
0, then, it can be shown that I> 2 = 2/a. Substituting and taking the 
hyperbolic tangent of both sides of equation (4.28), the result becomes 

F 3 =%tanh[-|=<-l)' e (l-e'^ aT1 ): . (4.29) 

Substituting this expression into equation (4.26) and simplifying, the 
resulting equation becomes 

F = sech[ ^ - ~ - (1-e 2ari )H , (4.30) 

1 2a 

Substituting equations (4.29) and (4.30) into equations (4. 21. a) and 
(4.21.b) s and substituting n - et and to = ae 9 the resulting closed-form 
solution for wave amplitudes and are 

c = e^hsechE ^ 1 ^ <l-e~ %a)t )]} (4. 31. a) 

1 2co 

— t / \ f 1 ~ ■ 

C, = { tanh[ e( — ( l-e~^ 1 * >t ) , (4.31.b) 

3 2 2(o 

To find expressions for and (Jig, substitute the relation that 
into equations (4.19.b) and (4.19.d) and integrate and evaluate the con- 
stants of integration with the initial conditions} the results are 
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h = ^30 = 2 ^10 _£7T C4,32) 

where 1 is a constant and $ is Zt\ radians out of phase with 2<J> . It 
10 _L 

can be seen that a special set of initial conditions is necessary to be 
consistent with this solution. A representative set is = ® 

which corresponds to Z - 0. 

Inspection of equations (4.31) reveals that the magnitude of 
continually decreases with time while the magnitude of first increases 
and then decreases. An interesting special case of equations (4.31) 
occurs in the absence of steady-state combustion (w = 0). The results of 
this case are 


C 1 = sech[ ( ] 

C 3 = | t anh [ ( ~ 0 — — ] ^ (4.33) 

These results show that a disturbance in the form of the first mode is 
transferred to the second mode as time increases. It is thought that this 
indicates the beginning of the steepening that leads to the formation of a 
shock wave. It can be seen that the presence of damping, in the form of 
steady-state combustion, inhibits this process. 

The second case to be investigated is that of standing waves with 
gas-dynamic nonlinearities neglected. To simulate the standing wave 
effect, let the amplitudes C 2 and C 4 and the phase angles <p 2 and 4> 4 equal 
zero. Again, this automatically satisfies four of the eight equations of 
(4.18). To achieve omission of gas-dynamic nonlinearities, let i = 0. 
Also, let the correction variable, K, be equal to zero for simplicity. 
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In doing so, the resulting equations, based upon equation C4.18), 
become 


dT + 1 + 2 ™ CC l C 3 sin(24, l ' ^3 ):1 = ° 

^ 4>i _ 

IT + 2 nu,[c 3 °° s(2 *l - V 1 = ° 
dC. i i 

d^r + 2 _aC 3 - 16 n " w[ " C l 2 Sin(2 ^l " *3 ):1 = 0 
d*, 1 C,* 

dT + il na,c "cJ cos(24, i " V 3 = 0 • 

The initial conditions imposed for this case are 


(4.34) 


c x (o) = 1 
c 3 (o) = 0 

*1<°> * ho 

*3 (0) * *30 . <4 ‘ 35> 

Let 2<j> - = (2 £ + 1) tt/2 , £ = 0, 1, 2 . . . . This implies that sin^ 

_ a) = (-1)^ and cos(2<}i 1 - <J> 3 ) = 0. Substituting into (4.34) and solving 
in the manner indicated previously one obtains expressions for the ampli- 
tudes for and which are 

_ -%U)t r -%0)t n -i ^ 

C = e tsec[— ne(-l) Cl-e )Aj 


(4. 36. a) 
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C„ = 


e'^ttanC^ 


(4. 36. b) 


3 " 2/2 


♦l = +10 


(4 . 36 . c) 


^3 “ 2<t> 10 ^ 2 17 


(4.36.d) 


where <j) in is constant and cj> is (2l+l)v/2 radians out of phase with 2$.. 

10 o 

As in the previous solution, special initial conditions are required to 
produce this solution, A representative set is $ = 0, <i> 30 = -v/2 9 which 

corresponds to t = 0. 

The secant and tangent both become infinite when their arguments 

take on the value ±ir/ 2. In (4. 36. a, b), the arguments of these functions 

3/2 

start at zero at t = 0 and have a maximum absolute value at ne/2 

Thus, If ne/2 3/ ^ < tt/ 2, the tangent and secant never become infinite and 

C and eventually decay to zero due to the influence of the exponential 
1 3 

3 /2 

function. This is a stable situation. If, on the other hand, ne/2 > 
tt/2, the tangent and secant become infinite at t^ = (2/w) |£n[l-2 2 Tr/(ne)] | 
causing C, and to become infinite. This is an unstable situation. 

1 o 

Thus, the boundary between stable and unstable behavior is indicated by 
the equation 


n e/2 3/2 = tt/2 . (4.37) 

The stability equation in the n-e plane has the form 

n = 2 K/e = 4.442/e . (4.38) 


This has the form of a rectangular hyperbola and is independent of o>. 
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For the case of traveling wayes, it is more convenient to work 
with the general perturbation equations expressed in modal amplitudes 
in terms of the real time variables, equation C4.13). To simulate the 
effect of spinning or traveling waves, let the following modal ampli- 
tudes be equal. These relations are 



It can be shown that substituting the relations (4.39) into equation 
(4.10), expressing the results in terms of the real time variables, sub- 
stituting these expressions into equation (3.18), and using appropriate 
multiple-angle formulas leads to 

<j)(0,t) = A^cos(t-O) - A2sin(t-0) + A^cos 2(t-0) 

-A^sin 2(t-0) + . . . . (4.40) 

which has the form of a sum of traveling waves. Substituting the expres- 
sions in (4.39) into equations (4.13), these eight equations reduce to 
four pairs of identical equations. The four independent equations listed 
below are 


dA i i 1 

— = e[- j aA 1 +^<A 2 -i(A 1 A 3 +A 2 A 4 )-na.(A 1 A 4 -A 2 A 3 ): 


dt 
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^ A 2 1 - 1 

— = e[- 2 aA 2 ~ 2 KA 1 - i (A 1 \- A 2 A 3) +nu)(A 2 A 4 +A l A 3 ):i 


3 1 - 1 1 

dt“ = eU " 2 aA 3 +4KA 4 + 4i( A 1 2 - A 2 2 )^CAi A 2 )3 

dA 4 1 - 1 1 - 

= E ^ _ 2 ctA 4 -4KA 3 + 2 i ^ A i A 2^“ gnu(A 1 2 -A 2 2 )] . (4.41) 


By making the substitution, we have reduced to a system of four equations 
and four unknowns. By solving for the modal amplitudes A , the modal 
amplitudes B_. are readily computed by using the relations of (4.39) to 
determine the entire nature of the wave form. 

For the case of traveling waves omitting gas-dynamic nonlinearities, 
let the amplitudes and A^ equal zero. Then set i, the gas-dynamic 
index, equal to zero. Again, for simplicity, let the correction variable, 
K, controlling physical chamber configurations, be zero. In doing so, in 
terms of the transformation variable, n, the resulting equations become 



nm[A 2 A 4 ] 


0 


1 - 1 - 

ITT + 2 a A 4 - 8 A 2 = ° 


(4.42) 


which is a system of two equations and two unknown modal amplitudes. To 
find an exact closed-form solution to these equations, let 


A 2 . 


\ ■ ^2 
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dA 

dn 


2 _-%an / 1 - 


= e 


(- 2 °^ F i + e 


-%an dF l 


dn 


kon, 1 - Nr , . „“%an 2 


2 r ■ t 5)f 2 * e 


dn 


(4.43) 


Using these transformations, the procedure for solution is exactly the 
same as for the standing wave case for both no combustion and no gas 
dynamics. The initial conditions for this case are 


a 2 (o) = 1 


A 4 (0) = o 


(4.44) 


Substituting the expressions of (4.42) into (4.41), the resulting equa- 
tions are 


s r - e- %on = o 





(4.45) 


with initial conditions 


F (0) = 1 


F 2 (0) = o 

Solving these equations in the manner outlined in the standing wave solu 
tions, the results are 
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F = S ec[44<l-e- %3n )] 

a 

F = (l-e' %5n )J . (4.46) 

2 2/2 2 5 

Expressing the results of (4.45) in terms of modal amplitudes by substi- 
tuting into (4.42), the resulting equations become 

A 2 = e-^seclff 2§ d-e- 43 ")] 

o 

A = L^l tan[~— (l-e" %OT1 )] . (4.47) 

4 2/2 2 a 


The results for traveling waves (4,47) are quite similar to the results 
for standing waves (4,36) for the case of no gas-dynamic nonlinearities. 
The same behavior can be expected as was discussed in the standing wave 
case about the nature of oscillation of the modal amplitudes. The only 
significant difference is the value to determine the boundary of stability 
for the interaction index governing the combustion terms. The stability 
condition for traveling waves is 

/2 TT 

— nE - — ^ (4.48) 


Thus, the equation of the stability boundary in the n-e plane is 


IT 



2.22 

e 


(4.49) 


Comparing equation (4.49) to (4.38) shows that the stability boundary for 
the interaction index is half as great for the traveling wave case as for 



61 

the standing wave case for any e. This will be yerified in a later pre- 
sentation of results of various numerical cases. 

For the case of traveling waves with no combustion, let the ampli- 
tudes and equal to zero. Then set n, the interaction index, equal 
to zero, and, again, let the correction variable K equal to zero. Sub- 
stituting into equations (4.40) and transforming into variable n> the 
results are 



A 1 A 3 


0 


dA 3 i - 1 2 

d n 2 ° A 3 " 4 A 1 


(4.50) 


with initial conditions 


A x (0) = 1 

a 3 (o) = 0 

which again is a system of two equations and two unknown modal amplitudes. 
To find an exact closed-form solution to these equations, use similar 
transformations as shown in (4.42). In doing so, and simplifying, the 
results are 


dF, 

— — + e 2<7r| F F 
dn 12 


o 


dF 


2 


dn 


I 2 = 

4 6 U 


0 


(4.51) 
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with initial conditions 

F 1 (0) = 1 
F 2 (0) = o . 

Solving these equations in the same manner as before, the results are 

F^ = sech[l/o(l-e ^ ar| )] 

F 2 = i tanh[l/a(l-e" %CTT1 )] . (4.52) 

Again, expressing the results of (4.51) in terms of the modal amplitudes 
of the form of equation (4.43), the resulting equations become 

= e ^ ari sechCl/a(l-e 2CTr| )H 

-%an i - 

Ag = — - — tanh[l/a(l-e 2an )] . (4.53) 

The results for the traveling waves (4.52) are similar to the results 
for standing waves (4.31) for the case of no combustion. A disturbance 
initially having the form of the first mode eventually is transformed into 
one having the form of the second mode. To compare these results for 
standing waves and traveling waves to the general perturbation equations, 
two computer programs were written (Appendices D and E) which numerically 
evaluate the modal amplitudes of various conditions for standing and 


traveling waves. 
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One last special case is an investigation of the effect of the 
correction variable K. In the special cases previously discussed, the 
correction variable K was set equal to zero. But, in this discussion, 
the correction variable K will be of primary importance in the equations. 
To start this analysis, refer to equations (3.21). Based upon these 
equations, impose the following conditions. First, neglect combustion 
effects (i.e., n = 0). Then, let us consider only the case of standing 
waves (i.e., g = g 2 = 0). Finally, let us neglect the steady state 
burning rate (i.e., o = 0) and assume that the terms multiplied by eK 
are larger than those multiplied by e above. This can be accomplished 
by writing 

K 1 = eK (4.54) 

and treating as a quantity of 0(1). Imposing the above conditions 
and substituting equation (4.54) into the equations (3.21), the result- 
ing equations become 

d 2 f, df df 

Cl+ V + f i + 2eCf 2 It" + f i dF- ] = 0 (4 ' 55 ' a) 

d 2f _ df 

[1+4K ] + 4f 0 - ef. -rr = 0 (4.55.b) 

1 dt 2 21 dt 

with initial conditions 


fjCo) = i 
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df. 


dt 


:o) = o 


f 2 (0) = 0 


df, 

J. 

dt" 


Ko) = o . 


First, assume a straightforward perturbation solution similar to the 
equations (4.6) except the functions are dependent upon the real time t. 
Substituting these assumed solutions into the equations and initial con- 
ditions of (4.55) and keeping terms of 0(1) and 0(e), the separated 
equations become 


d 2 f 


±° + 1 


dt 


2 (1+K^ 10 


f. = 0 


(4. 56. a) 


d 2 f 


20 4 


dt 


2 (1+4K 1 ) 20 


= 0 


(4. 56.b) 


d 2 f 


— + 


2 f df2 ° 


dt 2 


1+K. 11 ( 1+K ) L_ 20 dt 10 dt 


-] 


(4. 56. c) 


d 2 f. 




1 


9 ‘1+4K, 21 (1+4K ) 10 dt 

dt 2 1 1 


(4. 56.d) 


with initial conditions 


f io (0) = 1 


f u (0) = 0 
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df 


10 
dt 


< 0 ) = 0 


f 20 (0) 


= 0 


df 20 (0) 

dt 


= 0 


df 


11 


dt 


< 0 ) = 0 


f 21 (0) = 0 


df 2 i(0) 

dt 


= 0 


The first-order equations (4. 56. a and b) can be solved by assuming the 
usual assumed solution for linear differential equations. Doing this 
and applying the appropriate initial conditions, the results for the 
first-order terms are 


f 


10 


cos 


1+K, 


t 


f 


20 


0 


(4.57) 


Substituting (4.57) into the right-hand side of (4.56.c) the equation 
becomes a homogeneous linear differential equation. Solving in the 
usual manner and applying the appropriate initial conditions 

f =o (4.58) 

11 * 

Substituting (4.57) into the right-hand side of equation (4.56.d), the 
resulting equation becomes a linear differential equation with a particu- 
lar solution. By assuming an appropriate homogeneous and particular 
solution and evaluating the constants using the appropriate initial condi- 


tions , the result becomes 
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-A+4K 


'21 24K, 


1 2 I y— ■ 

— sin— t + A+K 


A+4K 


24K X 1 


sm 


VltK, 


t . (4.59) 


Therefore, substituting equations (4.57), (4.58), and (4.59) into the 
assumed perturbation solution and letting = eK, the resulting equations 
become 


f 


1 


cos 


1 

/ 1+Ke 


t 


(4. 60. a) 


f„ = 


-/I+Ke r /I+4KT . 


2 24K L / 1+Ke 


sm 


t - sin t] f (4*60. b) 

/l+4ek A+Ke 


Recall that in the two-variable perturbation method, and f 2 expressed 
in terms of the perturbation variables were 


f^ = A^(n) cos £ + B^(n) sin E, 


(4. 61. a) 


f 2 = Ag(n) cos 2E, + Bg(n) sin 2£ . (4.61.b) 

By transforming equation (4. 60. a) into perturbation variables and expand- 
ing the argument of the cosine function by the Taylor series and using 
appropriate sum and difference trigonometric identities f^ can be 
expressed as 


f = cos ^ Kri cos E, + sin — Kn sin £ 


(4.62) 


Therefore, comparing this to equation (4. 61. a), the functions and B 1 


must be 
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A 1 (n) 

B^n) 

By similar procedure, it can 
and comparing it to equation 

A 3 Cn) 

B 3 (n) 

To show the validity of equations (4.63) and (4.64), the problem is now 
solved using equations (4.12) which are derived from equations (.3.21) by 
the use of the two- variable perturbation method. To reproduce the condi- 
tions imposed on the problem just discussed, let there be no combustion 
(i.e., n = 0), let there be no steady -state burning rate (i.e., a = 0), 
and let there be only standing waves existing (i.e., = A^ = = B^ = 

0). Imposing these conditions on equations (4.12), the resulting equa- 
tions become 


= cos — Kn 


• 1 V 

sm j Kn . 


(4.63) 


be shown that evaluating equation (4.60.b) 
( 4. 61.b) , the results are 


24K 


'sin 4Kn - sin Kn3 


24K 


'cos Kn - cos 4KnU . 


(4.64) 


dA 


1 + j KB i + #: a i a 3 + B i B 3 i = 0 


dn 


dB i i i 

3T - 7 “l + f A 1 B 3 - Vs 3 ■ 0 


dA 


3 + 4KB 3 + ^CB^ - A^ 2 ] = 0 


dn 
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dB 3 1 

d^T - 4KA 3 - 4 A 1 B 1 = 0 • 


(4.65) 


In the previous solution it was assumed that the frequency correction 
terms were larger than the gas-dynamic nonlinearities. To be consistent 
with this assumption the following procedure is used. By a change of 
variable n = ?/K, equation (4.65) can be rewritten as 


dA l 1 1 

dc + 2 B 1 + 2K *- A l A 3 + B 1 B 3^ = 0 


dB 


1 1 Al + ^7 [A,B 0 - B. = 0 


d? 2 1 2K L T3 13- 


dA c 

dT 


f 4B 3 + IK CB 1 2 - V 3 = 0 


dB 3 1 

3T ' 4A 3 - 4K A 1 B 1 = 0 ’ 


(4.66) 


Assuming a straightforward expansion of the form 


A 1 A 10 + K A 11 + • • ' 


B 1 = B 10 + K B 11 + * ’ ’ 


A, = 


A 30 + K A 31 + 


B 30 + K B 31 + 


(4.67) 
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then substituting these expressions into the equations (_4. 
terms of 0(1) and 0(1/K), the resulting separated equatior 


dA 


10 4b,. = 0 


d£ 2 10 


dB 


10 4a. = 0 


d? 2 10 


^30 

“dT + 4B 30 ■ 0 


dB 30 „ 

— - 4A 30 = 0 


dA 


11 ,1b. 


d; 2 11 


" 2^ A 10 A 30 + B 10 B 30^ 


dB 


11 1 


d? 2 11 


" ^- A 10 B 30 " B 10 A 30 ] 


+ 4B 31 = - ^ B 10 Z - A 10 2] 


dB 


j~„ - 4^ - if a B 1 

dt 4A 31 4> A 1(T10 J 


with the initial conditions 


V 0) = 1 B io (0) = 0 


a i;l (o) = 0 


b 1 ]L Co) = o 


66) and keeping 
s become 

(4. 68. a) 
(4. 68 ,b) 
(4.68.c) 
(4.68.d) 
(4.68.e) 
( 4.68.f ) 
(4.68.g) 
(4.68.h) 
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A 30 (0) = ° B 30 (0) = ° 

A 31 (0) = 0 B 31 (0) = 0 * 


Since the first-order equations are coupled, differentiate equations 
(4. 68. a and c) once with respect to £ then substitute equations (4.68.b 
and d) into these equations resulting in 


d 2 A 




10 

2 


+ 


I A 

4 A 10 


0 


d 2 A 




30 

2 


16A 30 = ° 


(4.69) 


As can be seen, equations (4,69) are linear differential equations 
which can be evaluated by the usual manner. In doing so and applying 
the appropriate initial conditions, the resulting first-order modal ampli- 
tudes are 


. 1 1 „ 

A io = cos = cos 2^ 

fl 30 * 0 • ( “- 70) 

Knowing values for A and A , substitute these values into equations 

XU oU 

(4.68.b and d) and apply appropriate initial conditions. The results 
become 

• 1 . K 

B 1q = sin ^ = sin -^Kn 


B 


30 


= 0 . 


(4.71) 
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Substituting the results of (4.70) and (4.71) into the right-hand side of 
equations (4.68.e-h), the resulting equations become 


dA 


11 


d? 


-W 


(4. 72. a) 


dB ll i 
d c 2 A 11 


(4.72.b) 


dA 


31 


d£ 


+ 4B 


31 


1 

8 COS 




(4.72.c) 


dB 31 

"dT ' 4A 31 


1 . 

8 Sin 




(4.72. d) 


Since equations (4.72.C and d) are coupled, differentiate both equations 
once with respect to £ and substituting equations (4.72.C and d) into the 
appropriate terms of the new set of equations, the resulting equations are 


d 2 A 


31 


+ 16A 


d^ 


31 


sm £ 


d 2 B 


31 


+ 16B 




31 


5 

8 C ° S 


C 


(4.73) 


Equations (4.73) are a set of linear differential equations with homo- 
geneous and particular solutions. Solving these equations in the usual 
manner and using the appropriate initial conditions, the resulting modal 
amplitudes are 


31 24 


(sin4£ - sine) 



sin4Kq - sinKq) 
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31 24 


cos£ 


cos4?) = -^{cosKn - cos4Kn). 


(4.74) 


In a similar manner, the results for the modal amplitudes and can 
be determined to be 


B 1X = 0 (4.75) 

evaluated with the appropriate initial conditions. Therefore, substitut- 
ing the results of (4.70), (4.71), (4.73), and (4.74) into the assumed 
perturbation solution of (4.67), the resulting modal amplitudes become 

A 1 = cos y Kn + . . . 

= sin y Kn + . . . 


A3 = y— (sin 4Kn - sin Kn) + . . . 

B = -ri-(sin 4Kn - sin Kn ) + .... (4.76) 

o z4i\ 

It can be seen that equations (4.76) are identical to equations (4.63) and 
(4.64). This indicates that the two-variable method produces the correct 
solution. Equations (4.60) indicate that the presence of K changes the 
frequency of each of the first two acoustic modes and further renders the 
ratio of the second frequency to the first a non-integer number in general. 
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Equations (4.76) show how this effect manifests itself in the two-variable 
perturbation solution. 

These results can be used in another way. If the nonlinear terms 
are neglected in (4. 55. a), the results are 

d^f 

( 1+K ) =■ + f = 0 

1 dt 2 1 


(1+4K ) 



+ 4f , 


ef 


5 

1 dt 


0 


fjio) = 1 , 


df 1 (0) 

dt 


0 , 


f 2 (0) = 0, 


df 2 (0) 

dt 


0. (4.77) 


It can be easily shown that equations (4.60) constitute the exact solution 
of equation (4.77). If the corresponding terms are neglected in equations 
(4.65), the results are 


dA 


1 + 4 KB. =0 


d n 2 1 


dB 


1 - 4 KA, =0 


d n 2 1 


dA, 

dt 


+ 4KB 3 + |(B 1 2 - A x 2 ) = 0 


dB 3 1 

df ” 4KA 3 “ 4 A 1 B 1 = ° 


(4.78) 


where 
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AjtO) = 1 
B 1 (0) = 0 
AgCO) = 0 
Bg(O) = 0 • 

It can be shown that equations (4.76) are the exact solution of equation 
(4.78). These facts were used to check the accuracy of the computer pro- 
grams to be discussed later. 

In the remainder of this thesis, a comparison of the magnitudes 
of the modal amplitudes will be represented in graphical and tabular 
form. Under a given set of conditions, the acoustic modal amplitude pro- 
gram, the general perturbation program, and the analytical cases that 
were programmed will be used and results compared. Varying certain con- 
ditions will show their effect on the changes in magnitude of the modal 
amplitudes through a set range of time which is related to maintaining 
stability. From these various cases, it will be determined which param- 
eters and conditions have the greatest effect In changing modal ampli- 
tudes and which in turn affect the stability criteria for combustion 
by the methods discussed above. 


Chapter 5 


DISCUSSION AND PRESENTATION OF RESULTS 

In this chapter , results are presented both in graphical and 
tabular form which are representative of the results generated by the 
programs listed In the Appendices B through E. From these representative 
sets of results, basic observations will be made to observe which 
parameters or conditions have the greatest effects on the problems of 
stability . 

In Figures 3 and 4, modal amplitudes F^ and F^ are graphically 
represented versus time for a stable standing wave case. For these 
figures, F^O) = 0, F 1 '(0) = 1, F 2 (0) = 0, F 2 '(0) = 0, 6 1 (0) = 0, G^CO) = 
0, G 2 (0) = 0, G 2 ’(0) = 0, n = 35, i = 1, K = 0, e = 0.1 and u> = 0.1. 

The step size used was 0.1. Experimentation showed that this was a small 
enough step size to produce accurate results and was used throughout. 

From these figures, one notices that both the first and second order modal 
amplitudes decrease In amplitude with increasing time. Also, T^ 9 the 
second order modal amplitude, tends to oscillate at twice the frequency of 
F^. These figures are based upon one set of parametric values; however, 
these figures represent qualitatively the results obtained using a wide 
variety of initial conditions and parametric values. In Figures 5 through 
8, modal amplitudes F^, F 2 , , and G 2 are graphically represented versus 

time for a stable traveling wave case. For these figures, F^(0) = 0, 
F 1 , (0) = -1, F 2 (0) = 0, G 1 (0) = 1, G 1 , (0) = 0, G 2 (0) = 0, G 2 '(0) = 0, 
n = 15, i = 1, K = 0, to = 0.1 and e = 0.1. The general shape of the 
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curves and the relative frequencies of oscillation are qualitatively 
similar to the stable standing wave case. 

In Figures 9 and 10, modal amplitudes F 1 and F 2 are graphically 
represented versus time for an unstable standing wave case with the 
same conditions as the stable case except that n = 50. As can be seen, 
the maximum amplitude of F^ starts to decrease then increase dramatically 
for increasing time. The maximum amplitude of F^ increases continuously. 

In Figures 11 through 14, modal amplitudes F^, F 2 , , and G 2 are 

represented versus time for an unstable traveling wave case. Again, the 
conditions are the same as for the stable traveling wave case except that 
n = 30. Drastic increases in amplitudes are observed for all the modal 
amplitudes shown as time increases. The behavior is similar to the 
unstable standing wave case. The period of time for traveling waves to 
become unstable is about one-half the period of time for standing waves 
to become unstable. Thus, it seems that traveling waves are less 
stable than are standing waves. 

In Tables 1 and 2 , a comparison of results is presented for modal 
amplitudes F^ and F 2 for a stable standing wave case. For these cases, 
F 1 (0) = 0, F 1 '(0) = 1, F 2 (0) = 0, F 2 '(0) = 0, G a (0) = 0, G ± ' ( 0 ) = 0, 

G 2 (0) = 0, G 2 ’(0) = 0, n = 60, e = 0.1, and m = 0.1. These tables 
quantitatively show the effect of neglecting gas dynamic non-linearities 
on the accuracy of the computations. Also, a comparison can be made 
between the exact solution method (Appendix B program) and the perturbation 
solution method (Appendix C program). From Table 1, one can observe that 
the effect of neglecting gas-dynamic nonlinearities is small where 
quantitatively comparing values of the modal amplitude F^. Even though, 
quantitatively, the values for the exact solutions and perturbation 
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Table 1. Comparison of Results for F^ Showing Effects of Gas Dynamic 
Index (i) - (Fj_ = 0, F^ = 1, F 2 = 0, F 2 ' = 0, G 1 = 0, = 0 , G 2 = 0, 

G 2 T = 0) - Stable Cases (n = 60) - Standing Waves 



i 

K 

= 1 
= 1 

i = 0 
K = 1 

t 

Exact 

Perturbation 

Exact 

Perturbation 


Solution 

Solution 

Solution 

Solution 

0.2 

0.19699 

0.18712 

0.19699 

0.18702 

0.4 

0.38335 

0.36426 

0.38336 

0.36386 

0.6 

0.55252 I 

0.52540 

0.55259 

0.52462 

0.8 

0.69856 

0.66518 

0.69885 

0.66400 

1.0 

0.81627 

0.77905 

0.81719 

0.77758 

1.2 

0.90132 

0.86340 

0.90354 

0.86186 

1.4 

0.95043 

0.91572 

0.95489 

0.91443 

1.6 

0.96159 

0.93461 

0.96936 

0.93399 

1.8 

0.93432 

0.91986 

0.94632 

0.92040 

2.0 

0.86985 

0.87244 

0.88656 

0.87466 

2.2 

0.77125 

0.79443 

0.79242 

0.79884 

2.4 

0.64330 

0.68895 

0.66783 

0.69602 

2.6 

0.49211 

0.56003 

0.51822 

0.57016 

2.8 

0.32469 

0.41247 

0.35021 

0.42593 

3.0 

0.14827 

0.25167 

0.17115 

0.26856 

3.2 

-0.03017 

0.08340 

-0.01142 

0.10366 

3.4 

-0.20430 

-0.08633 

-0.19032 

-0.06300 

3.6 

-0.36859 

-0.25158 

-0.35912 

-0.22566 

3.8 

-0.51829 

-0.40659 

-0.51242 

-0.37879 

4.0 

-0.64917 

-0.54604 

-0.64583 

-0.51727 

4.2 

-0.75738 

-0.66519 

-0.75583 

-0.63655 

4.4 

-0.83921 

-0.76006 

-0.83953 

-0.73278 

4.6 

-0.89118 

-0.82756 

-0.89450 

-0.80297 

4.8 

-0.91029 

-0.86556 

-0.91868 

-0.84504 

5.0 

-0.89444 

-0.87297 

-0.91049 

-0.85789 

5.2 

-0.84300 

-0.84979 

-0.86909 

-0.84143 

5.4 

-0.75726 

-0.79704 

-0.79483 

-0.79656 

5.6 

-0.64071 

-0.71679 

-0.68959 

-0.72514 

5.8 

-0.49885 

-0.61200 

-0.55708 

-0.62988 

6.0 

-0.33868 

-0.48646 

-0.40279 

-0.51427 

6.2 

-0.16794 

-0.34466 

-0.23365 

-0.38244 

6.4 

0.00574 

-0.19159 

-0.05742 

-0.23901 

6.6 

0.17554 

-0.03259 

0.11809 

-0.08892 

6.8 

0.33582 

0.12685 

0.28578 

0.06271 

7.0 

0.48221 

0.28126 

0.43974 

0.21080 

7.2 

0.61134 

0.42538 

0.57552 

0.35042 

7.4 

0.72042 

0.55435 

0.68999 

0.47703 

7.6 

0.80682 

0.66387 

0.78099 

0.58653 

7.8 

0.86778 

0.75032 

0.84693 

0.67550 

8.0 

0.90042 

0.81090 

0.88644 

0.74121 
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Table 2. Comparison of Results for F 2 Showing Effects of Gas Dynamic 
Index (i) - (Fj_ = 0, = 1, ? 2 = 0, F 2 ' = 0, G x = 0, G ’ = 0, G 2 = 0, 

G 2 ' = 0) - Stable Cases (n = 60) - Standing Waves 



i = 1 
K = 1 

i = 0 
K = 1 

t 

Exact 

Solution 

Perturbation 

Solution 

Exact 

Solution 

Perturbation 

Solution 

0.2 

0.00012 

-0.00485 

0.00003 

-0.00253 

0.4 

0.00113 

-0.01309 

0.00043 

-0.00938 

0.6 

0.00422 

-0.02223 

0.00205 

-0.01865 

0.8 

0.01060 

-0.02944 

0.00602 

-0.02784 

1.0 

0.02110 

-0.03215 

0.01336 

-0.03424 

1.2 

0.03582 

-0.02854 

0.02471 

-0.03553 

1.4 

0.05397 

-0.01795 

0.03497 

-0.03023 

1.6 

0.07375 

-0.00104 

0.05816 

-0.01798 

1.8 

0.09260 

0.02022 

0.07742 

0.00026 

2.0 

0.10749 

0.04283 

0.09521 

0.02232 

2.2 

0.11543 

0.06317 

0.10863 

0.04514 

2.4 

0.11407 

0.07764 

0.11493 

0.06516 

2.6 

0.10215 

0.08319 

0.11198 

0.07891 

2.8 

0.07988 

0.07788 

0.09878 

0.08357 

0 

00 

0.04909 

0.06128 

0.07573 

0.07745 

C\l 

00 

0.01309 

0.03461 

0.04478 

0.06031 

3.4 

-0.02375 

0.00071 

0.00929 

0.03352 

3.6 

-0.05653 

-0.03632 

-0.02639 

-0.00006 

00 

CO 

-0.08051 

-0.07164 

-0.05749 

-0.03640 

4.0 

-0.09180 

-0.10031 

-0.07945 

-0.07083 

4.2 

-0.08798 

-0.11801 

-0.08865 

-0.09863 

4.4 

-0.06850 

-0.12165 

-0.08296 

-0.11571 

4.6 

-0.03490 

-0.10993 

-0.06211 

-0.11921 

4.8 

0.00924 

-0.08353 

-0.02786 

-0.10793 

5.0 

0.05868 

-0.04515 

0.01610 

-0.08257 

5.2 

0.10711 

0.00077 

0.06458 

-0.04571 

5.4 

0.14796 

0.04864 

0.11144 

-0.00152 

5.6 

0.17534 

0.09236 

0.15041 

0.04468 

5.8 

0.18493 

0.12615 

0.17593 

0.08715 

6.0 

0.17466 

0.14534 

0.18393 

0.12042 

6.2 

0.14513 

0.14699 

0.17255 

0.14004 

6.4 

0.09957 

0.13035 

0.14243 

0.14314 

6.6 

0.04352 

0.09700 

0.09680 

0.12888 

6.8 

-0.01592 

0.05073 

0.0411 

0.09857 

7.0 

-0.07104 

-0.00295 

-0.01771 

0.05557 

7.2 

-0.11448 

-0.05745 

-0.07205 

0.00488 

7.4 

-0.14027 

-0.10594 

-0.11471 

-0.04743 

7.6 

-0.14458 

-0.14223 

-0.13988 

-0.09498 

7.8 

-0.12628 

-0.16158 

-0.14384 

-0.13189 

8.0 

-0.08716 

-0.16128 

-0.12556 

-0.15352 
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solutions are not exactly the same, the order of magnitude and behavior 
of results is similar. From Table 2, the same observations can be made 
for the behavior of F 2 . There is, however, more error, quantitatively, 
between the results for exact and perturbation methods and a region of 
qualitative inaccuracy between the exact and perturbation solutions exists 
near t = 0. This takes the form of a difference in sign of F 2 between 
results from the exact solution as compared to the perturbations solution. 
This discrepency occurred also in the other calculations performed (not 
shown) and will be discussed in more detail later in this chapter. 

3h Tables 3 and 4, a comparison of results is presented for modal 
amplitudes F^ and F 2 for a stable standing wave case. The initial 
conditions for the results in these tables are F^(0) = 0, F^ (0) - 1, 

F 2 (0) = 0, F 2 f (0) = 0, G 1 (0) = 0, G 1 T (0) = 0, G 2 (0) = 0, G 2 '(0) = 0, 
n = 40 , e =0.1, and u> = 0.1. However, these tables quantitatively 
present the effect of deviations of the ratio of the second acoustic 
frequency to the first from the integer value of 2 (this is controlled 
by the parameter K) . These results show that solutions for finite values 
of K are qualitatively similar to those for K = 0. This indicates that 
the ratio of the second acoustic frequency to the first does not have 
to be an integer in order to produce the type of behavior observed here . 

A ratio near an integer value will lead to similar results. Tables 3 
and 4 also allow a comparison to the results generated by the program 
in Appendix D for the approximate analytical solution (4.31). These 
results presented in the last column of Tables 3 and 4 can be compared 
to the fourth column in each of these tables to determine the accuracy 
of (4.31). These comparisons present further evidence that the neglect 
of gas dynamic nonlinearities does not have an important qualitative 


effect . 
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Table 3. Comparison of Results for Showing Effects of the Correction 
Variable (K) - (F^ = 0, F^' = 1, F2 = 0, F2' = 0 , G-j_ = 0, G^' = 0, 

G2 = 0, G 2 ' = 0) - Stable Case (n = 40) - Standing Waves 



i = 
K = 

1 

1 

n 11 

1 

0 

■nQI 
1 1 

t 

Exact 

Perturbation 

Exact 

Perturbation 



Solution 

Solution 

Solution 

Solution 

1 

0.2 

0.19699 

0.18707 

0.19670 

0.19678 

0.19671 

0.4 

0.38335 

0.36396 

0.38172 

0.38210 

0.38186 

0.6 

0.55254 

0.52460 

0.54795 

0.54890 

0.54843 

0.8 

0.69867 

0 . 66358 

0.68903 

0.69093 

0.69029 

1.0 

0.81667 

0.77635 

0.79957 

0.60302 

0.80234 

1.2 

0.90247 

0.85935 

0.87537 

0.88124 

0.88075 

1.4 

0.95312 

0.91014 

0.91361 

0.92306 

0.92304 

1.6 

0.96699 

0.92744 

0.91310 

0.92740 

0.92820 

M- 

00 

0.94390 

0.91120 

0.87443 

0.89468 

0.89666 

2.0 

0.88523 

0.86255 

0.80001 

0.82676 

0.83027 

2.2 

0.79388 

0.78374 

0.69393 

0.72690 

0.73221 

2.4 

0.67411 

0.67805 

0.56163 

0.59954 

0.60682 

2.6 

0.53129 

0.54969 

0.40948 

0.45016 

0.45947 

00 

CM 

0.37154 

0.40358 

0.24432 

0.28905 

0.29627 

3.0 

0.20133 

0.24520 

0.07300 

0.11102 

0.12386 

3.2 

0.02712 

0.08039 

- 0.09782 

- 0.06484 

- 0.05085 

3.4 

- 0.14491 

- 0.08439 

- 0.26190 

- 0.23547 

- 0.22097 

3.6 

- 0.30901 

- 0.24474 

- 0.41335 

- 0.39410 

- 0.37989 

00 

00 

- 0.45992 

- 0.39354 

- 0.54662 

- 0.53453 

- 0.52152 

4.0 

- 0.59287 

- 0.52617 

- 0.65660 

- 0.65139 

- 0.64055 

4.2 

- 0.70357 

- 0.63813 

- 0.73872 

- 0.74028 

- 0.73261 

4.4 

- 0.78821 

- 0.72574 

- 0.78928 

- 0.79800 

- 0.79446 

4.6 

- 0.84362 

- 0.78623 

- 0.80583 

- 0.82264 

- 0.82407 

4.8 

- 0.86746 

- 0.81785 

- 0.78760 

- 0.81364 

- 0.82069 

5.0 

- 0.85849 

- 0.81988 

- 0.73571 

- 0.77179 

- 0.78490 

5.2 

- 0.81682 

- 0.79266 

- 0.65317 

- 0.69920 

- 0.71852 

5.4 

- 0.74409 

- 0.73757 

- 0.54464 

- 0.59917 

- 0.62456 

5.6 

- 0.64352 

- 0.65697 

- 0.41585 

- 0.47611 

- 0.50704 

5.8 

- 0.51967 

- 0.55405 

- 0.27311 

- 0.33525 

- 0.37090 

6.0 

- 0.37816 

- 0.43279 

- 0.12272 

- 0.18253 

- 0.22172 

6.2 

- 0.22516 

- 0.29772 

0.02937 

- 0.02427 

- 0.06554 

6.4 

- 0.06694 

- 0.15382 

0.17767 

0.13305 

0.09140 

6.6 

0.09056 

- 0.0063 

0.31714 

0.28307 

0.24291 

6.8 

0.24192 

0.13957 

0.44298 

0.41977 

0.38308 

7.0 

0.38242 

0.27867 

0.55061 

0.53776 

0.50650 

7.2 

0.50795 

0.40617 

0.63563 

0.63247 

0.60850 

7.4 

0.61493 

0.51774 

0.69408 

0.70029 

0.68527 

7.6 

0.70017 

0.60967 

0.72786 

0.73878 

0.73407 

7.8 

0.76090 

0.67898 

0.72020 

0.74667 

0.75327 

8.0 

0.79476 

0.72355 

0.68609 

0.72400 

0.74234 
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Table 4. Comparison of Results for F 2 Showing Effects of the Correction 
Variable (K) - (F^ = 0, ' = 1, F 2 = 0, F 2 = = = 

G 2 = 0, G 2 ' = 0) - Stable Case (n = 40) - Standing Waves 



i = 
K = 

1 

1 

i = 1 
K = 0 

^ H- 
II II 
O O 

t 

Exact 

Perturbation 

Exact 

Perturbation 

Analytic 


Solution 

Solution 

Solution 

Solution 

Solution 

0.2 

0.00011 

-0.00400 

0.00015 

-0.00419 

-0.00192 

0.4 

0.00099 

-0.00996 

0.00136 

-0.01035 

-0.00696 

0.6 

0.00354 

-0.01599 

0.00479 

-0.01597 

-0.01337 

0.8 

0.00860 

-0.02013 

0.01141 

-0.01856 

-0.01885 

1.0 

0.01665 

-0.02071 

0.02149 

-0.01626 

-0.02113 

1.2 

0.02761 

-0.01669 

0.03443 

-0.00838 

-0.01856 

1.4 

0.04072 

-0.00793 

0.04866 

0.00438 

-0.01059 

1.6 

0.05458 

0.00479 

0.06180 

0.01995 

0.00208 

1.8 

0.06726 

0.01983 

0.07113 

0.03529 

0.01749 

2.0 

0.07663 

0.03491 

0.07416 

0.04696 

0.03278 

2.2 

0.08072 

0.04750 

0.06910 

0.05185 

0.04472 

2.4 

0.07803 

0.05521 

0.05549 

0.04793 

0.05038 

2.6 

0.06793 

0.05622 

0.03439 

0.03468 

0.04776 

2.8 

0.05081 

0.04957 

0.00841 

0.01343 

0.03626 

0 

00 

0.02816 

0.03542 

-0.01867 

-0.01282 

0.01697 

3.2 

0.00244 

0.01507 

-0.04246 

-0.03978 

-0.00745 

3.4 

-0.02323 

-0.00917 

-0.05871 

-0.06265 

-0.03313 

3.6 

-0.04544 

-0.03425 

-0.06410 

-0.07695 

-0.05561 

3.8 

-0.06099 

-0.05679 

-0.05687 

-0.07938 

-0.07067 

4.0 

-0.06735 

-0.07356 

-0.03727 

-0.06854 

-0.07507 

4.2 

-0.06307 

-0.08196 

-0.00764 

-0.04530 

-0.06722 

4.4 

-0.04802 

-0.08040 

0.02778 

-0.01278 

-0.04759 

4.6 

-0.02355 

-0.06854 

0.06351 

0.02408 

-0.01872 

4.8 

0.00767 

-0.04745 

0.09363 

0.05926 

0.01509 

5.0 

0.04187 

-0.01945 

0.11279 

0.08668 

0.04845 

5.2 

0.07470 

0.01210 

0.11719 

0.10127 

0.07573 

5.4 

0.10173 

0.04330 

0.10529 

0.09987 

0.09206 

5.6 

0.11915 

0.07013 

0.07824 

0.08191 

0.09414 

5.8 

0.12425 

0.08905 

0.03976 

0.04959 

0.08095 

6.0 

0.11586 

0.09745 

-0.00434 

0.00764 

0.05395 

6.2 

0.09462 

0.09402 

-0.04711 

-0.03736 

0.01700 

6.4 

0.06286 

0.07894 

-0.08154 

-0.07807 

-0.02424 

6.6 

0.02438 

0.05389 

-0.10178 

-0.10757 

-0.06318 

6.8 

-0.01630 

0.02182 

-0.10409 

-0.12058 

-0.09336 

7.0 

-0.05322 

-0.01335 

-0.08754 

-0.11433 

-0.10952 

7.2 

-0.08237 

-0.04728 

-0.05423 

-0.08921 

-0.10854 

7.4 

-0.09956 

-0.07576 

-0.00908 

-0.04879 

-0.09001 

7.6 

-0.10237 

-0.09523 

0.04092 

0.00093 

-0.05639 

7.8 

-0.09018 

-0.10322 

0.08777 

0.05190 

-0.01267 

8.0 

-0.06428 

-0.09872 

0.12375 

0.09586 

0.03434 
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In Tables 5 and 6, a comparison of results are presented for modal 
amplitudes and for an unstable standing wave showing the effect of 
neglecting gas-dynamic nonlinearities. It can be seen that the gas 
dynamic nonlinearities have little qualitative effect on the results. 

In Tables 7 and 8, a comparison of results are presented for modal 
amplitudes F^ and F^ for an unstable standing wave case showing the effects 
of K. The results for zero and non-zero are qualitatively similar. 

These tables are representative of the cases that were investigated 
in the course of this research. Only cases involving standing waves were 
presented. The same behavior, however, can be observed for the cases 
involving traveling waves. 

In Table 9, a comparison of stability boundaries is presented 
based upon the interaction index (n) which is a measure of the strength 
of the combustion process. For standing waves and the given conditions 
shown, the stability limit for a process with gas dynamic nonlinearities 
considered and K = 0 is between 45-50. When both gas dynamic non- 
linearities and the correction variable are considered, the stability 
limit is increased to 67.5-69. Finally, when considering only the 
correction variable with no gas-dynamic non-linearity effect, the stability 
limit is 72-72.5. The results show that the neglect of gas dynamic 
non linearities slightly underestimates the stability boundary and that 
the increasing K increases the stability limit. 

In Table 10, a comparison of stability boundaries is presented 
based upon the interaction index for traveling waves. These results provide 
additional confirmation of the conclusions discussed in the previous 
paragraph and also illustrate the fact that standing waves are roughly 
twice as stable as traveling waves. This is consistent with the 



Table 6. Comparison of Results for F 2 Showing Effect of the Gas 
Dynamic Index (i)- (F^=0,F^' =1 9 F 2 =0, F 2 ' =0, G^=0, 

Gi' = 0, G 2 = 0, G 2 ' = 0) - Unstable Case - (n = 75) - Standing Waves 



i = 1 
K = 1 

i = 0 
K = 1 

t 

Exact 

Perturbation 

Exact 

Perturbation 


Solution 

Solution 

Solution 

Solution 

0.2 

0.00013 

-0.00548 

0.00003 

-0.00316 

0.4 

0.00124 

-0.01544 

0.00054 

-0.01172 

0.6 

0.00473 

-0.02691 

0.00257 

-0.02333 

0.8 

0.01210 

-0.03643 

0.00752 

-0.03484 

1.0 

0.02443 

-0.04077 

0.01670 

-0.04288 

1.2 

0.04197 

-0.03749 

0.03089 

-0.04728 

1.4 

0.06387 

-0.025518 

0.04994 

-0.03793 

1.6 

0.08806 

-0.005425 

0.07263 

-0.02258 

1.8 

0.11143 

0.02061 

0.0966 

0.00033 

2.0 

0.13027 

0.049007 

0.1186 

0.02813 

2.2 

0.14082 

0.07537 

0.13504 

0.05696 

2.4 

0.14005 

0.09513 

0.14236 

0.08237 

2.6 

0.12628 

0.10426 

0.13796 

0.09995 

2.8 

0.09967 

0.10004 

0.12064 

0.10606 

3.0 

0.06242 

0.08152 

0.09109 

0.09851 

3.2 

0.01860 

0.04984 

0.05192 

0.07688 

3.4 

-0.02636 

0.008166 

0.0076 

0.04283 

3.6 

-0 . 06636 

-0.03858 

-0.03628 

-0.00096 

3.8 

-0.09548 

-0.08436 

-0.07362 

-0.046816 

4.0 

-0.10889 

-0.12285 

-0.09877 

-0.09134 

4.2 

-0.10354 

-0.14826 

-0.10733 

-0.12758 

4.4 

-0.07873 

-0.15626 

-0.09688 

-0.15016 

4.6 

-0.03636 

-0.14461 

-0.06742 

-0.15524 

4.8 

0.01909 

-0.11361 

-0.02153 

-0.14104 

5.0 

0.08104 

-0.06617 

0.03583 

-0.10827 

5.2 

0.14148 

-0.00758 

0.09778 

-0.06010 

5.4 

0.19206 

0.055118 

0.15638 

-0.00186 

5.6 

0.22528 

0.11399 

0.20358 

0.05954 

5.8 

0.23565 

0.16121 

0.23235 

0.11651 

6.0 

0.22067 

0.19016 

0.23779 

0.16168 

6.2 

0.18130 

0.19633 

0.21796 

0.18888 

6.4 

0.12201 

0.17793 

0.17434 

0.19397 

6.6 

0.05019 

0.13633 

0.11182 

0.17547 

6.8 

-0.02482 

0.07595 

0.03808 

0.13477 

7.0 

-0.09302 

-0.02526 

-0.03737 

0.07614 

7.2 

-0.14503 

-0.07144 

-0.10448 

0.006209 

7.4 

-0.17338 

-0.140235 

-0.15407 

-0.06680 

7.6 

-0.17347 

-0.19368 

-0.17900 

-0.13401 

7.8 

-0.14428 

-0.22450 

-0.17519 

-0.18700 

8.0 

-0.08859 

-0.22808 

-0.14216 

-0.21888 



Table 5. Comparison of Results for F a Showing Effect of the Gas Dynamic 
Index (i) - (Fj = 0, F x ' = 1, F 2 = 0, F 2 ' = 0, G ± = 0. G ± ' = 0 , G 2 = 0, 
Go' = O') - Unstable Case (n = 75) - Standing Waves 




Exact 

Solution 


0.19699 

0.38335 

0.55250 

0.69847 

0.81592 

0.90030 

0.94803 

0.95671 

0.92557 

0.85572 

0.75039 

0.61484 

0.45593 

0.28147 

0.09945 

-0.08272 

-0.25864 

-0.42307 

-0.57177 

-0.70107 

-0.80745 

-0.88721 

-0.93641 

-0.95124 

-0.92862 

-0.86710 

-0.76771 

-0.63434 

-0.47368 

-0.29439 

-0.10598 

0.08258 

0.26374 

0.43198 

0.58355 

0.71613 

0.82779 

0.91635 

0.97878 

1.01111 


Perturbation 

Solution 


0.18717 

0.38215 

0.52613 

0.69244 

0.80625 

0.86736 

0,92126 

0.94184 

0.92873 

0.88276 

0.80579 

0.70079 

0.57161 

0.42287 

0.25987 

0.08829 

-0.08588 

-0.25669 

-0.41827 

-0.56512 

-0.69228 

-0.79547 

-0.87123 

-0.91703 

-0.93138 

-0.91381 

-0.86492 

-0.78635 

-0.68071 

-0.55149 

-0.40299 

-0.22558 

-0.06815 

0.10715 

0.28009 

0.44478 

0.59627 

0.72909 

0.83897 

0.92224 


Exact 

Solution 


0.19699 
0.38336 
0.55258 
0.69881 
0.81700 
0.90294 
0.95337 
0.96695 
0.94003 
0.87584 
0.77580 
0.64411 
0.48677 
0.31118 
0.12554 
-0.06194 
-0.24371 
-0.41342 
-0 . 56608 
-0.69790 
-0.80594 
-0.88767 
-0.94056 
-0.96203 
-0.94961 
-0.90146 
-0.81711 
-0.69815 
-0.54872 
-0.37551 
-0.18717 
0.00672 
0.19691 
0.37555 
0.53682 
0.67707 
0.79428 
0.88737 
0.95535 
0.99664 


Perturbation 

Solution 


0.18703 

0.36402 

0.52513 

0.66518 

0.77974 

0.86535 

0.91952 

0.94089 

0.92921 

0.88529 

0.81106 

0.70941 

0.58410 

0.43963 

0.28107 

0.11388 

-0.05620 

-0.22350 

-0.38241 

-0.52773 

-0.63473 

-0.75934 

-0.83825 

-0.88902 

-0.91020 

-0.90121 

-0.86255 

-0.79563 

-0.70274 

-0.58703 

-0.45231 

-0.30299 

-0.14391 

0.01979 

0.18290 

0.34025 

0.48686 

0.61813 

0.72996 

0.81891 
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Table 6. Comparison of Results for F 2 Showing Effect of the Gas 
Dynamic Index (i) - (F* = 0, Fp’ = 1 , F 2 = 0, F 2 ' = 0, G 1 = 0, 

G^' = 0, G 2 = 0, G 2 ' = 0) - Unstable Case - (n = 75) - Standing Waves 



i = 1 
K = 1 


O vH 
it it 

t 

Exact 

Solution 

Perturbation 

Solution 

Exact 

Solution 

Perturbation 

Solution 

0.2 

0.00013 

-0.00548 

0.00003 

-0.00316 

0.4 

0.00124 

-0.01544 

0.00054 

-0.01172 

0.6 

0.00473 

-0.02691 

0.00257 

-0.02333 

0.8 

0.01210 

-0.03643 

0.00752 

-0.03484 

1.0 

0.02443 

-0.04077 

0.01670 

-0.04288 

1.2 

0.04197 

-0.03749 

0.03089 

-0.04728 

1.4 

0.06387 

-0.025518 

0.04994 

-0.03793 

1.6 

0.08806 

-0.005425 

0.07263 

-0.02258 

1.8 

0.11143 

0.02061 

0.0966 

0.00033 

2.0 

0.13027 

0.049007 

0.1186 

0.02813 

2.2 

0.14082 

0.07537 

0.13504 

0.05696 

2,4 

0.14005 

0.09513 

0.14236 

0.08237 

2.6 

0.12628 

0.10426 

0.13796 

0.09995 

2. 8 

0.09967 

0.10004 

0.12064 

0.10606 

3.0 

0.06242 

0.08152 

0.09109 

0.09851 

3.2 

0.01860 

0.04984 

0.05192 

0.07688 

3.4 

-0.02636 

0 . 008166 

0.0076 

0.04283 

3 . 6 

-0.06636 

-0.03858 

-0.03628 

-0.00096 

3 . 8 

-0.09548 

-0.08436 

-0.07362 

-0.046816 

4.0 

-0.10889 

-0.12285 

-0.09877 

-0.09134 

4.2 

-0.10354 

-0.14826 

-0.10733 

-0.12758 

4.4 

-0.07873 

-0.15626 

-0.09688 

-0.15016 

4.6 

-0.03636 

-0.14461 

-0.06742 

-0.15524 

4.8 

0.01909 

-0.11361 

-0.02153 

-0 . 14104 

5.0 

0.08104 

-0.06617 

0.03583 

-0.10827 

5.2 

0.14148 

-0.00758 

0.09778 

-0.06010 

5.4 

0.19206 

0.055118 

0.15638 

-0.00186 

5.6 

0.22528 

0.11399 

0.20358 

0.05954 

5.8 

0.23565 

0.16121 

0.23235 

0.11651 

6.0 

0.22067 

0.19016 

0.23779 

0.16168 

6.2 

0.18130 

0.19633 

0.21796 

0.18888 

6.4 

0.12201 

0.17793 

0.17434 

0.19397 

6.6 

0.05019 

0.13633 

0.11182 

0.17547 

6. 8 

-0.02482 

0.07595 

0.03808 

0.13477 

7.0 

-0.09302 

-0.02526 

-0.03737 

0.07614 

7.2 

-0.14503 

-0.07144 

-0.10448 

0.006209 

7.4 

-0.17338 

-0.140235 

-0.15407 

-0.06680 

7.6 

-0.17347 

-0.19368 

-0.17900 

-0.13401 

7. 8 

-0.14428 

-0.22450 

-0.17519 

-0.18700 

8.0 

-0.08859 

-0.22808 

-0.14216 

-0.21888 



97 


Table 7. Comparison of Results for F a Showing the Effects of the 
Correction Variable (K)-(F^=0, F^ = 1. ^2 = ^2 - 0, G-^-0, 

G t ' = 0, G 2 = 0, G 2 ' = 0) - Unstable Cases (n = 70) - Standing Waves 



i = 3 
K = 1 


n n 

•H ^ 

1 

0 

i = 0 
K = 0 

t 

Exact 

Solution 

Perturbation 

Solution 

Exact 

Solution 

Perturbat ion 
Solution 

Analytic 

Solution 

0.2 

0.19699 

0.187155 

0.19670 

0.19687 

0.19675 

0.4 

0.38335 

0.36462 

0.38172 

0.38261 

0.38217 

0.6 

0.55251 

0.52587 

0.54791 

0.55028 

0.54942 

0.8 

0.69850 

0.66615 

0.68878 

0.69371 

0.69249 

1.0 

0.81604 

0.78072 

0.79865 

0.80767 

0.80630 

1.2 

0.90066 

0.86596 

0.87280 

0.88811 

0.88696 

1.4 

0.94887 

0.91928 

0.90775 

0.93226 

0.93184 

1.6 

0.95842 

0.93925 

0.90166 

0.93878 

0.93468 

1.8 

0.92863 

0.92554 

0.85477 

0.90772 

0.91061 

2.0 

0.86067 

0.87904 

0.76962 

0.84061 

0.84612 

2,2 

0.75771 

0.80166 

0.65109 

0.74802 

0.74903 

2.4 

0.62481 

0.69648 

0.50590 

0.61107 

0.62333 

2.6 

0.46860 

0.56735 

0.34194 

0.45807 

0.47407 

2.8 

0.29658 

0.41903 

0.16733 

0.26849 

0.30714 

3.0 

0.11649 

0.25679 

-0.01027 

0.10616 

0.12906 

3.2 

-0.06441 

0.08639 

-0.18416 

-0.05313 

-0.05327 

3.4 

-0.23973 

-0.08618 

-0.34860 

-0.25962 

-0.23280 

3.6 

-0.40414 

-0.25497 

-0.49854 

-0.42949 

-0.40262 

3.8 

-0.55318 

-0.41415 

-0.62912 

-0.58151 

-0.55621 

4.0 

-0.68301 

-0.55829 

-0.73533 

-0.70960 

-0.68766 

4.2 

-0.78998 

-0.68250 

-0.81197 

-0.80859 

-0.79192 

4.4 

-0.87041 

-0.78258 

-0.85399 

-0.87449 

-0.86494 

4.6 

-0.92055 

-0.85525 

-0.85736 

-0.90458 

-0.90388 

4.8 

-0.93689 

-0.89811 

-0.82008 

-0.89755 

-0.90718 

5.0 

-0.91671 

-0.90981 

-0.74306 

-0 . 85356 

-0.87463 

5.2 

-0.85886 

-0.89005 

-0.63045 

-0.80213 

-0.80736 

5.4 

-0.76447 

-0.83962 

-0.48925 

-0.66264 

-0.70785 

5.6 

-0.63724 

-0.76029 

-0.32818 

-0.52309 

-0.57982 


-0.48340 

-0 . 65486 

-0.15634 

-0.36111 

-0.42807 

1 E ! ' 

-0.31100 

-0.52691 

0.01809 

-0.18308 

-0.25835 


-0.12888 

-0.38083 

0.18852 

-0.07645 

-0.07712 

Kn 

0.05445 

-0.22156 

0.34995 

0.10764 

0.10864 


0.23170 

-0.05445 

0.49844 

0.37422 

0.29176 

■ 

0.39723 

0.11485 

0.63023 

0.46032 

0.46509 

7.0 

0.54707 

0.28072 

0.74089 

0.69003 

0.62177 

7.2 

0.67846 

0.43759 

0.82500 

0.81052 

0.75551 

7.4 

0.78917 

0.58027 

0.87632 

0.89873 

0.86085 

7.6 

0.87687 

0.70394 

0.88888 

0.95061 

0.93330 

7.8 

0.93868 

0.80457 

0.85847 

0.96357 

0.96962 

00 

o 

0.97109 

0.87877 

0.78419 

0.93647 

0.96786 
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Table 8. Comparison of Results for F 2 Showing the Effects of the 
Correction Variable (K) - (F^ = 0, F^ = 1, F 2 = 0, F 2 = 0, - 0, 

G 1 t = o, G 2 = 0, G 2 ’ = 0) - Unstable Cases (n = 70) - Standing Waves 



i = 0 
K = 0 


Exact 

Solution 


0.00013 
0.00120 
0.00456 
0.01160 
0.02332 
0.03992 
0.06057 
0.08330 
0.10517 
0.12271 
0,13242 
0.13150 
0.11839 
0.09328 
0.05822 
0.01701 
-0.02528 
-0.06295 
-0.09049 
-0.10337 
-0.09872 
-0.07590 
-0.03665 
0.01487 
0.07254 
0.12595 
0.17637 
0.20783 
0.21824 
0.20519 
0.16949 
0.11514 
0.04883 
-0.02090 
-0.08483 
-0.13429 
-0.16228 
-0.16442 
-0.13960 
-0.09020 


Perturbation 

Solution 


-0.00527 

-0.01466 

-0.02534 

-0.03411 

-0.03789 

-0.03451 

-0.02299 

-0.00396 

0.02047 

0.04692 

0.07126 

0.08922 

0.097153 

0.11651 

0.07467 

0.04468 

0.00567 

-0.03775 

-0.07994 

-0.11504 

-0.13779 

-0.1442 

-0.13261 

-0.1032 

-0.05894 

-0.00478 

0.05266 

0.10617 

0.14862 

0.17409 

9.17865 

0.16089 

0.12229 

0.06705 

0,00165 

-0.06594 

-0.12723 

-0.17433 

-0.20087 

-0.21857 


Exact 

Solution 


0.00017 

0.00165 

0.00619 

0.01544 

0.03024 

0.05013 

0.07306 

0.09555 

0.11324 

0.12170 

0.11751 

0.09917 

0.06770 

0.02679 

-0.01761 

-0.05827 

-0.08789 

-0.10041 

-0.09221 

-0.06290 

-0.01566 

0.04299 

0.10404 

0.15731 

0.19315 

0.20424 

0.18706 

0.14284 

0.07764 

0.00147 

-0.0733 

-0.13397 

-0.16972 

-0.17346 

-0.14308 

-0.08205 

0.00095 

0.09312 

0.17949 

0.24522 


Perturbation 

Solution 


-0.00562 

-0.01557 

-0.02603 

-0.03276 

-0.03218 

-0.02232 

-0.01862 

0.021911 

0.04915 

0.07264 

0.08675 

0.087046 

0.071414 

0.04063 

-0.001408 

-0.04828 

-0.09196 

-0.12426 

-0.13827 

-0.12986 

-0.09854 

-0.04785 

0.01493 

0.07988 

0.13602 

0.17311 

0.183575 

0.163922 

0.11562 

0.04509 

-0.03703 

-0.11757 

-0.182899 

- 0.22122 

-0.22468 

-0.19089 

-0.12367 

-0.03255 

0.06839 

0.16276 
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Table 9. Comparison of Stability Boundaries Based on the Interaction 
Index (n) - (F-^ = 0, = 1, F 2 = 0, F 2 = 0, G^ = 0, G-^ - 0, 

G 2 = 0, G 2 ' = 0) - Standing Waves - Epsilon - 0.1 


Gas Dynamic Index 
Correction Variable 

Stability 

Boundaries 

Exact Solution 
n - Stable - Unstable 

Perturbation Solution 
n - Stable - Unstable 

i = 1 

67.5 - 69 

67.5 - 69 

11 



i = 0 

72 - 72.5 

72.5 - 73 

K = 1 




i = 1 
K = 0 


45 - 50 


45 - 50 
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Table 10. Comparison of Stability Boundaries Based on the Interaction 
Index (n) - (F a = 0, T ± ' = -1, F 2 = 0, F 2 ' = 0, G 1 = 1, G ± ' =0, 

G 2 = 0, G 2 ' = 0) - Traveling Waves - Epsilon - 0.1 
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approximate analytical stability equations (4.31) and (4.52). The 
perturbation method tends to predict slightly higher stability limits 
than the exact solution method for both standing and traveling waves. 
Within the accuracy of the tabulated values, this is apparent only in 
the first two rows of Table 10 . 

In Table 11, a comparison of the effect of different initial 
conditions imposed on the stability boundaries for both standing and 
traveling waves is presented. From the results of two sets of initial 
conditions for each case, it can be seen that the varying of initial 
conditions has no significant effect on the stability boundaries for 
both standing waves or traveling waves. 

In Table 12 ? the variation of the stability limit with e is 
presented for standing waves. From Table 12, the results show that the 
smaller the term epsilon the greater the stability limit. Therefore, 
the order term has a significant effect on the interaction index. In 
Chapter 4, a relation was proposed for the case of i = 0 and K = 0 
which was n = C/e where C is a constant. Assuming the validity of the 
relation, the values for this constant are given for each given epsilon 
and interaction index. This shows that, in general, C is a weak function 
of e . 

In Table 13, a comparison of the effect of e is presented for 
traveling waves when both gas dynamic nonlinearities and correction 
variables are considered. Again, the results show that the smaller the 
term epsilon, the greater the stability limit. The perturbation method 
again predicts slightly greater stability limits than does the exact 
solution method. Therefore, again, the order term has a strong effect 
concerning the stability of combustion. 
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Table 11. Comparison of the Effect of Different Initial Conditions 
Imposed for Standing and Traveling Waves for i = 1 and K = 1 
Epsilon = 0.1 

(a) Standing Waves - 1. = 0, F^ T = 1, F 2 = 0, F 2 * = 0 

= 0 , G^ = 0 , G 2 = 0 5 G 2 = 0 

2. F 1 = 1, = 0, F 2 = 0, F 2 T = 0 

Gj_ = 0 5 — 0, G 2 ~ 0, G 2 ~ 0 


Initial Condition 
Sets 

Stability Boundaries 

Exact Solution 
n - Stable - Unstable 

Perturbation Solution 
n - Stable - Unstable 

1. 

67.5 - 69 

67.5 - 69 

2. 

65 - 70 

65 - 70 


(b) Traveling Waves - 1. F^ = 0, Fj/ = -1, F 2 = 0, ¥ 2 * - 0 

G^ - 1 5 G^ = 0 , G 2 ” 0 5 G 2 = 0 


2 . F 1 = 1 , F t ’ = 0 , F 2 = 0, F 2 ’ = 0 
G-^ — 0 ? G^ ~ —1 9 G 2 = 0 , G 2 — 0 


Initial Condition 
Sets 

Stability Boundaries 

Exact Solution 
n - Stable - Unstable 

Perturbation Solution 
n - Stable - Unstable 

1. 

27.5 - 28 

31.5 - 32 

2. 

27.5 - 28.5 

31 - 31.5 
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Table 12. Comparison of the Effects of the Order Term Epsilon - 


(F 1 = 0, F 1 ' = 1, F 2 = 0, F 2 ' = 0, Gj. = 0, Gi* = 0, G 2 = 0, 
G 2 * = 0) - Standing Waves - when i = 1, K = 1 


Epsilon 

Stability Boundaries 


Exact Solution 
n - Stable - Unstable 

Perturbation Solution 
n - Stable - Unstable 

Constant 
C = ns 

0.05 

107.5 - 110 

107.5 - 110 

5.5 

0.1 

67.5 - 69 

67.5 - 69 

6.9 

0.2 

48.5 - 49.5 

■ 

48.5 - 49.5 

9.8 




















Table 13. Comparison of the Effects of the Order Term Epsilon (F^ = 0, F^ ' = -1, F 2 = 0, 
F 2 1 = 0, = 1, 6 j' = 0, G 2 = 0, 62 ' =0) Traveling Waves - when i = 1, K = 1 


Epsilon 

Stability Boundaries 

Exact Solution 
n - Stable - Unstable 

Perturbation Solution 
n - Stable - Unstable 

Constant C = ne 
(Exact Solution) 

Constant C = ne 
(Perturbation Sol) 

0.05 

51 - 52 

52 - 53 

2.575 

2.625 

0.1 

27.5 - 28 

31.5 - 32 

2.775 

3.175 

0.2 

15.5 - 16 

18.5 - 19.5 

3.15 

3.75 
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Thus, from these representative tables of results , it is observed 
that the correction variable is important in the stability of standing 
waves, but does not play a major role in the stability of traveling waves. 
It is observed that the gas dynamic nonlinearities seem to have little 
influence on the stability of either standing or traveling waves. It 
is observed that initial conditions of the modal amplitudes have little 
or no influence in the stability of either standing or traveling waves. 

And finally, it is observed that the order term epsilon and* the inter- 
action index governing the strength of combustion in the process are 
strongly coupled thus affecting the limits of stability. 

Before completing this chapter, it is desired to investigate the 
sign discrepency mentioned previously between the exact and perturbation 
solutions for f which occur near t = 0. For simplicity, it will be 
assumed that i = K = 0 and that for t << 1 the first modal amplitude can 
be represented with sufficient accuracy by f^ = sint. Then, the 
equation for f^ will be solved and the result simplified for t << 1. 

This will be done first for oj = 0 and then for w 1 0. For oj = 0 , 

(3.21) leads to 


d 2 f 

Z_ -j- up = \ 

dt2 + Hr 2 


% eoon 


- cos2tj 


with initial conditions 


(5.1) 


f 2 (0) = 0 


f 2 '(0) = 0. 
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Evaluating the homogeneous and particular solutions by the usual manner 
and evaluating the constants, the results become 


f. = eton 1 - cos2t - t sin2t . 
2 16 L J 


(5.2) 


In terms of the perturbation parameters (4.1), equation (5.2) can be 
written as 


= 16 [' 


e(l - cos2£) - n sin2£ 


■> • 


(5.3) 


To the order of approximation £ which the perturbation solution should 
model, equation (5.3) becomes 


f 


2 


1 

16 


wnn sin2£ + 0(e), 


(5.4) 


By expanding equation (5.2) into a Taylor series expansion of three terms, 
equation (5.2) becomes 


f 


2 


1 

24 


eomt 


+ . 


(5.5) 


which is always positive. 

Therefore, the exact method for small time will yield f^ modal 
amplitude always as a positive quantity. 

By imposing identical conditions to the perturbation equations 
(4.12), the result becomes 


1 

16 
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dB 

2 . = 

dn 


un 


(5.6) 


with the condition 


B 2 (0) = 0. 

Solving equation (5.6), 

B 2 = - wnn . (5.7) 


Recalling that f 2 = B 2 sin2£, the result becomes 

f 0 = - — wnn sin2£ + 0(e) 

/ lb 


(5.8) 


which is identical to the result of equation (5.4) for the wave equation 
solution. Thus, the perturbation method gives the correct result. It 
can be seen that for t << 1 the exact solution predicts a positive f 2 
and by inspection of equation (5.8), the perturbation method predicts a 
negative . This is precisely the behavior observed in the numerical 
solutions. 

For u> t 0, a similar analysis can be performed. The appropriate 
equation for f 2 is now 


d 2f df r 

_2. t i t 4 f 2 = V»nl 1 


cos2t 


] 


(5.9) 
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with conditions 
f 2 (0) = 0 
f 2 T (0) = 0. 

Solving the homogeneous and particular solution by the usual manner and 
evaluating the appropriate constants the result becomes 


* -0)/2t f 1 - / A 6 - or ^ 

f 2 = e - is £wn cos l — 2 1 


en 

16 


' 8 - 


£5 2 


sm 


/l6~ 


0)^ 


t + 


U) Z . 


1 

16 


1 

£wn * 8 


en sin2t 




(5.10) 


Expanding (5.10) for small a) into the appropriate Taylor series, expanding 
and neglecting terms of 0(ca) leads to 

f 2 = 16 

which is identical to (5.2). 

By imposing the identical conditions on the perturbation equation 
(4.12), the resulting equation become 


cos2t - t sin2t 


(5.11) 


dB . 

dn 


+ ^ oB 2 = 


1 

16 


can 


(5.12) 


with the condition 
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B 2 (0) = 0. 

Solving equation (5.12) by the usual manner, evaluating the constants, 
and transforming the perturbation variables to real time variables 

f 2 = - 2£[i _ e -*] sin 2t. (5.13) 

This is always negative for t << 1, Expanding the exponential function 
by the Taylor series expansion and neglect terms of o(w) leads to 


f 2 = n sin 25 + 0(e) 


(5.14) 


which is identical to (5.8). 

To observe the behavior of equation (5.10) for small time, 
expand this equation into a Taylor series of 0(t 4 ). Expanding and 
grouping terms according to their order of magnitude, the terms of 
0(1), 0 ( t ) , 0(t 2 ), 0(t 3 ) vanish. Therefore, f^ is comprised of terms 
from 0(t 4 ) which is 


f 


2 


- 4 r -4 

enwt 3ca w 

24 L 1 + 8 64 , 


(5.15) 


Again, for any small time t, f^ is always positive since t is always 
positive. Neglecting higher powers of o>, the resulting equation becomes 
equation (5.5) for the undamped case. Again it can be seen that the 
exact and perturbation methods predict opposite signs for when t << 1. 
These results are based on approximations and cannot be considered 
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definitive. They do, however, lend plausibility to the numerical results 
discussed earlier. It is believed that this sign discrepancy is due 
to the inability of the perturbation solution to accurately represent 
the exact solution for t << 1 and not due to any error in the computer 
program used to compute the perturbation solution. 


Chapter 6 


CONCLUSION AND RECOMMENDATIONS 

The primary objective of this presentation has been the development 
of analytical techniques to solve the problem of combustion Instabilities 
occurring in an annular combustion chamber. The analytical techniques 
used were the modified Galerkin method applied to the acoustic wave 
equations which yielded a set of time-dependent modal amplitude equations 
and the two-variable perturbation method which yield a set of time- 
dependent equations which approximated the behavior of the first set of 
equations. Both methods produced results which were relatively easy to 
apply and used the Runge-Kutta algorithm which required little computation 
time. An alternative approach to solve this problem would be a finite 
difference approach. However, difficulties can be foreseen In the 
development of the finite difference equations modelling the problem 
along with the complications occurring due to the boundary conditions of 
the problem. Thus, the benefits of the methods discussed in this thesis 
can be appreciated. 

From the numerical and graphical presentation of results In Chapter 
5, the following observations can be made. First, the effect of the gas- 
dynamic nonlinearities seems to be small in both methods of analysis for 
velocity sensitive combustion. This point can be observed from a 
quantitative comparison of the tabular results or by observing the effects 
of this condition on the stability boundaries. Second, the effect of the 
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correction variable modelling the physical boundaries of the chamber seems 
to have a significant effect in both methods of analysis for velocity 
sensitive combustion. By including the effect of this correction variable, 
a significant increase occurs in the interaction index which is the 
criteria for the stability of the system. However, this effect seems to 
be more significant for the standing wave case than the traveling wave 
cases. The effects of initial conditions for the time dependent equations, 
the numerical value for the burning rate and step size of integration, 
seem to have very little significance in the measure of the stability 
limits of velocity sensitive combustion. However, the order term epsilon 
has a strong effect upon the stability of the problem. This is to be 
excepted since the order term is the measure of the effect of non- 
linearities occurring in the system. The increase in this value corresponds 
to a decrease in the stability limit which is physically reasonable. 

In this study, the effect of time delay of the combustion process 
was neglected. However, time delay has been found in other studies to 
be an important phenomena in correctly modelling the actual problems of 
velocity sensitive combustion. It is recommended that this effect can 
be incorporated by Including the corresponding terms with j = 1 in the 
acoustic wave equations (3.20). A corresponding set of perturbations can 
then be derived to account for time delay and both these equations and 
equations (3.20) can be numerically evaluated by modifing the existing 
Runge-Kutta programs presented in the Appendices. It is also recommended 
that an experimental program be developed to measure the effects of 
velocity sensitive combustion in an annular combustion chamber. Once 
achieving this goal, one could correlate the measurement results to the 
analytical results that have been presented to ascertain the validity of 
this analysis. 
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Since instability of combustion is sensitive to small changes in 
engine geometry and operating conditions, a particular engine must be 
subjected to a large number of firings before its designers can say 
confidently that it is free from instability. With a large engine such 
testing can account for a substantial part of development costs. Herein 
lies the importance of devising reliable theories of instability and 
inexpensive tests of a propellant's acoustical characteristics. Until 
instability of combustion is understood well enough so that it can be 
eliminated while an engine is in the design stage, rocket engines must 
continue to be intensively tested for stability — particularly when 
the lives of astronauts will eventually depend on safe, reliable 
operation of the engine [17]. 
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GENERAL TIME DELAY FUNCTION 


The development and nature of the time-delay function is of the 
same form of the convolution integral for impulse response in vibration 
theory. The general form of the time delay function is 


w(t) 


t 

| j(t - o 
0 



(A.l) 


A simple illustration of the time delay function is in the case of a 

finite step function J(t). 

J(t) 

1 1 


T 

(some specific time constant) 
Figure Al. Step Function J(t) 


t 


From the figure, the step function J(t) is defined as 


J(t) 


1 t < T 
] 0 t > T 


(A. 2) 


Therefore, substituting some time delay (t - 5) for time t, the result is 
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Figure A2. Step Time Delay Function J(t - £) 


Substituting into the general time-delay integral the particular step 
function in terms of the non-dimensional variable £ 


ft - t du> 


0 ) 


(t) = 


ft du. 


0 dC . 


d? 


dK 


t - T 


(A 


Therefore, simplifying equation (A. 3) 


w(t) = w Q (t) - w Q (t “ T ) 


(A 


where w^Ct) is a generalized function of time and w Q (t - t) is 


functional time delay. 




APPENDIX B 

RUNGE-KUTTA PROGRAM OF THE MODAL 
AMPLITUDE WAVE EQUATIONS 
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APPENDIX D 


PROGRAM FOR EXACT SOLUTION OF 


STANDING WAVE CASE 
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ACOUSTIC PRESSURE DERIVATION 


To calculate expressions for acoustic pressure, recall equation 
(2.48) which stated 


P = P 




(F.l) 


This equation represents the unsteady state deviations of acoustic 
pressure. When expanding equation (F.l) into a Taylor series expansion, 
the resulting equation becomes 


P = P = 1 




+ . 


(F.2) 


Recall that the steady state solution was represented in equation (2.35) by 


P = e 




When expanding (F.3) into its Taylor series expansion, the result becomes 


P = 1 


W2 



2 


+ . . 


(F . 4 ) 
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where p is the steady state acoustic pressure. Therefore, the difference 
in general acoustic pressure and steady state pressure can be expressed 
by subtracting equation (F.4) from (F.2). For this investigation, a 
restriction on the velocity potential <J> was that it was a function of 0 
and t only. In doing this, the pressure difference equation becomes 



Using the same Fourier series expansion for the velocity potential <p as 
expressed in equation (3.18), the acoustic pressure difference equation 
(F.5) can be expressed in terms of the product of modal amplitudes and 
trignometric function in the transverse 0 direction. Substituting the 
appropriate forms of equation (3.18) into equation (F.5) and simplyfying, 
the resulting pressure difference equation become 


£ - £ = 
e 


df , r I 

dt^ + E P f l f 2 + B l g 2> + 15 l 


dg dg df df 
dt dt + dt dt 


COS0 


df 

2 

dt 


+ £ 


2 2 (V lf f fis V 

-*<«! - y -tar 1 ) 


cos20 


dg 

dt 


+ e 


" (f l g 2 " f 2 g 


1 } + h [< 


'df dg dg df 

dt dt dt dt J 


sinG 


dg 

2 

dt 


+ e 


/df \ /dg 


1 '2 f 1 g 1 + h \ dt i \ dt 



sin20 , 


(F. 6 ) 
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Since the coefficients in equation (F.6) are functions of time only, 
these coefficients have been included in the calculations of the program 
in Appendix B. Thus, for any given angle 0, values for the modal 
amplitude at any given time range can be calculated therefore determining 
the acoustic pressure difference of that desired location. 



